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I . INTRODUCTION 


1.0 Background 

To date, successful satellite capture and recovery/repair 
operations have been executed by Space Shuttle astronauts using the 
Remote Manipulating Arm and the Manned Maneuvering Unit. 1 * 2 However, 
both of these valuable tools are very much range-limited. Furthermore, 
attempts to capture uncontrolled satellites within the range of the 
Orbiter can lead to the dangerous exposure of the Space Shuttle crew 
members and/or the Orbiter to a possibly violently gyrating, 
uncontrolled satellite. For both of these reasons most future capture 
and recovery missions will likely be the domain of remotely piloted 
Orbital Maneuvering Vehicles (OMVs). 3 * 4 

The potential usefulness of OMVs to affect the capture and 
retrieval of satellites should be significantly increased by adaptable, 
yet accurate, simulations of the target satellite's attitude motions. 
Two primary mathematical models are needed for such simulation 
capability. First, a long-term model is needed for predicting a 
satellite's attitude motion as much as six months or more in advance. 
Second, a "short-term" model adaptable and rigorous enough for use in 
simulators to train the OMV "pilots" for specific satellite recoveries, 
is needed. 

Although the problem of modeling a spacecraft's motion about its 
center of mass is one which arose because of space flight, the study of 


1 


2 


spacecraft attitude dynamics is basically a contemporary application of 
classical theoretical mechanics. The analytical approach to solving the 
differential equations of rotational motion of a satellite is typified 
by Jacobi's classical solution for the free rotational motion of a 
triaxial rigid body in terms of elliptic integrals and functions. 
Although analytical solutions have been obtained for only relatively 
simple models of satellites, the more recent wide availability of 
digital computers has made it possible to use straightforward numerical 
integration of differential equations of motion for models composed of 
rigid and flexible bodies, subjected to various torques and forces. 

1.1 The Capture and Retrieval Problem 

A time period of several months may elapse between the failure of a 
satellite (for example, a control system failure, power system failure, 
maneuvering propellant depletion, or other non-catastrophic type 
problems which might precipitate the need for capture and retrieval) and 
the first possible opportunity, or "window," to capture and retrieve it. 
Hence, an apparently necessary capability is that of accurately 
predicting the evolution of the attitude motion of a disabled satellite 
several months into the future, for the purposes of capture feasibility 
analyses and mission planning. In addition to long-term prediction 
capability, short-term attitude motion simulation capability is needed 
to train remote operators. Since all satellites which may be the 
objects of capture and retrieval missions will most certainly not fit 
into a single category of satellite or attitude motion, simulation 
capability must include more than just a single-axis-of-rotation , 
single-rigid-body model. The satellite model should be rigorous enough 


3 


to include the presence of spinning rotors or reaction wheels 
distributed symmetrically, or asymmetrically . The flexibility 
of the satellite might also be included. 

Another problem involves the quantity and quality of data on the 
states of the satellite's orbital and attitude motions. These, of 
course, are functions of the satellite's "mode of failure." There are 
three general modes, or scenarios, of "failure" which may precipitate a 
recovery. The first scenario involves a satellite which has not lost 
its telemetry, but is uncontrollable. In this case, ground controllers 
would, through the satellite's telemetry data, have at least partial 
attitude motion information. The second scenario involves a satellite 
which is both uncontrolled and has lost its telemetry capability. Such 
a satellite's state was possibly known prior to failure, but additional 
information can only be obtained through remote observation, if the 
satellite is observable. The third scenario involves space debris which 
are, of course, neither controllable nor capable of telemetry. Thus, 
the debris' rotational states are unknown except for estimates based on 
any observations which might have been possible. Therefore, in the 
capture of debris, any of a wide variety of rotational states may be 
encountered. 5 

1.2 Scope of this Investigation 

The primary purpose of this research is to provide mathematical 
models which may be used in investigation of various aspects of the 
remote capture and retrieval of uncontrolled satellites . Emphasis has 
been placed on analytical models; however, to verify analytical 
solutions, numerical integration must be used. Also, for satellites of 
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of certain types, numerical integration may be the only practical or 
perhaps the only possible method of solution. First, to provide a basis 
for analytical and numerical work, uncontrolled satellites were 
categorized using criteria based on (1) orbital motions, (2) external 
angular momenta, (3) internal angular momenta, (4) physical 
characteristics and (5) the stability of their equilibrium states. 
Chapter II deals specifically with categorization. 

Second, several analytical solutions for the attitude motions of 
satellite models were compiled, checked, corrected in some minor 
respects and their short-term prediction capabilities were investigated. 
These models should be useful in studying the attitude behavior of 
satellites which have considerable angular momentum, and in "driving" 
short-term simulation, during which human operators try to capture 
satellites. Single-rigid-body, dual-spin and multi-rotor configurations 
are treated. The models are described in Chapter III and Appendices A, 

B and C of this report . Copies of computer codes for evaluating the 
solutions are being supplied separately. 

Third, as indicated above, to verify the analytical models and to 
see how the "true" motion of a satellite which is acted upon by 
environmental torques differs from its corresponding torque-free motion, 
a numerical simulation code was developed. This code contains a 
relatively general satellite model and models for gravity-gradient and 
aerodynamic torques. The spacecraft physical model for the code and the 
equations of motion are given in Chapter IV. The two environmental 
torque models are described later in Chapter V and Appendix D. 
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Fourth, the use of torque-free analytical solutions to represent 
satellite attitude motion is considered in Chapter VI. Analytical 
results and numerical results, including gravity-gradient and 
aerodynamic torque effects for satellites are presented. 

Fifth, there are cases in which analytical solutions are not known 
for the unperturbed attitude motion of a satellite, but available 
integrals of the motion and extensive numerical results indicate that 
the motion will in a great many cases be almost periodic. To model the 
motion in such cases, a semi-analytic method was developed. This method 
is based on concepts from the Generalized Method of Averaging (GMA) , but 
its application is primarily numerical. The method is described in 
Chapter VII. It has been used to predict satellite attitude motion over 
long periods of time. Results obtained using the semi-analytic method 
and direct numerical integration are presented in Chapter VII. 

Finally, Chapter VIII contains conclusions and recommendations. 

It should be noted here that this effort did not include an 
investigation of the effects of internal dissipation of energy on 
satellite motion. This was considered to be beyond the scope of this 
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II. CATEGORIZATION OF SATELLITES 


2.0 Rationale 

Artificial earth satellites can be categorized according to (1) the 
orbital motion of their centers of mass, (2) their gross (external) 
rotational motion, (3) their internal motion, (4) their individual 
physical characteristics and (5) the stability of their equilibrium 
states. 5 The motion of a satellite's center of mass about the earth is 
of importance in determining whether the satellite may be reached with 
the means available. If the satellite is accessible, the orbital 
characteristics are needed in modeling the effects of environmental 
torques. The initial rotational motion of the satellite about its 
center of mass as well as its internal angular momentum determine the 
satellite's initial rotational kinetic energy and total rotational 
angular momentum. Using these and satellite physical data, one can 
predict the general principal characteristics of the uncontrolled 
attitude motion which may evolve. Examples of this uncontrolled 
attitude motion include continuous, steady rotation and librating 
motion. 

2.1 Categories 

The block diagram shown in Fig. 1 depicts categories of 
uncontrolled satellites based upon orbital criteria. 5 The first two 
categories, LEO (Low Earth Orbit) and HEO (High Earth Orbit), are based 
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Fig. 1. Block Diagram of Categories of Uncontrolled Satellites. 
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upon the perigee altitude of the satellite's orbit. A LEO is considered 
one which has a perigee altitude of less than 600 km. For those periods 
of time during which the satellite operates in the LEO altitude range, 
it is subjected to significant aerodynamic torques as well as other 
environmental torques. For that portion of the satellite's orbit in 
which it moves in the HEO altitude range , aerodynamic torques become 
much less significant than gravity-gradient torques. 

The next two categories, NCO (Near Circular Orbit) and E0 (Elliptic 
Orbit), determine the satellite's orbital motion with respect to 
altitude and time. An elliptic orbit of significant eccentricity 
subjects the satellite to substantial periodic changes in the 
environmental torques which affect the satellite's attitude motion. On 
the other hand, for a satellite moving in a circular, or near circular, 
orbit, the magnitudes of the environmental torques to which it is 
subjected remain in a narrow range. 

Uncontrolled satellites also can be placed into several 
subcategories based on their rotational motion about their respective 
centers of mass. These subcategories include High Kinetic Energy (HKE) , 
Low Kinetic Energy (LKE) , High Angular Momentum . (HAM) , and Low Angular 
Momentum (LAM) . The HKE satellites possess rotational kinetic energy 
which is orders of magnitude larger than the maximum possible potential 
energy due to gravity-gradient torque acting upon the satellite. This 
high kinetic energy state could be due to a high "gross" rate of 
rotation; that is, a high rotational rate of the entire body, and/or due 
to large amounts of internal angular momentum and internal rotational 


kinetic energy. Satellites in the HAM subcategory possess rotational 
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angular momentum, _H, due either to gross or internal relative motion of 
sufficient magnitude that the ratio of the magnitude of the applied 
torque to the magnitude of H is small « Obviously, the subcategories HAM 
and HKE are not mutually exclusive. In point of fact, they overlap a 
great deal. 

Other subcategories are defined to deal more specifically with the 
inertia characteristics of the individual satellites. Two of these 
subcategories include Asymmetric Inertia (AI) ellipsoid, and Symmetric 
Inertia (SI) ellipsoid. Still other subcategories may be based upon the 
initial orientation of the satellite’s body-fixed axes with respect to 
the radius vector between the center of the earth and the satellite's 
center of mass. This orientation is of central importance in the 
determination of the torque on the satellite which in turn are modeled 
in stability analyses. The stability of its initial motion has a direct 
bearing upon the later motion of a satellites under the effects of 
aerodynamic and/or gravity-gradient torques. Those subcategories, 
associated with aerodynamic torques include AS, ANS, AUS (Aerodynamical- 
ly Stable, Neutrally Stable, and Unstable initial states, respectively). 
Similarly, for Gravity Gradient stability, we have the subcategories 
GGS , GGNS, and GGUS . 5 In all these the type of stability considered is 
static stability rather than dynamic stability. 

Finally, satellites can be categorized according to their construc- 
tion and their components. Such categories relate directly to the 
choice of the physical model and will approximate a satellite best. 
Examples are dual-spin satellites and satellites which contain reaction 


wheels for attitude control. 


III. ANALYTICAL MODELS FOR TORQUE-FREE 


ROTATIONAL MOTION OF SATELLITES 


3.0 Introduction 

The kinds of satellites to be considered in this investigation as 
typical candidates for capture and retrieval include single-rigid-body 
and dual-spin satellites as well as satellites which contain multiple 
spinning rotors (reaction wheels). In this Chapter three torque-free 
analytical models for these three types of satellites are discussed. 

Such models serve at least three purposes. First, if its rotational 
angular momentum is not too small, over short time periods such as that 
required for capture, the actual attitude motion of a satellite is 
essentially the same as its torque-free motion. Second, to predict 
satellite motion over longer periods of time, perturbation theories for 
satellite attitude motion, 6-10 analogous to those for satellite orbital 
motion may be developed on the basis of analytical solutions for 
torque-free motion. Third, analytical solutions provide a means for 
checking the accuracy of numerical integration procedures. 

3.1 Single-Rigid-Body Satellite Model 

As the name Implies , the single-rigid-body model consists of a 
single, rigid, asymmetric body. Using this single-rigid-body model, and 
general perturbation methods. Beletski, 9 Crenshaw and Fitzpatrick, 6 
Cochran , 1 and Liu and Fitzpatrick, 10 have developed theories to predict 
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the attitude motion In the presence of gravity-gradient torques. If 
the satellite modeled (see Fig. 2) is rotating rapidly about its center 
of mass, then its short-term motion is essentially classical "free- 
Eulerian" motion. This "free-Eulerian" motion may be described by 
an exact analytical solution which involves elliptic functions and 
integrals. 11 * 12 Appendix A contains the details of this mathematical 
model . 

3.2 Dual-Spin Satellite Model 

Satellites which contain a single spinning rotor and a nominally 
despun platform, or "carrier” body, are called "dual-spin” satellites 
(see Fig. 3). More specifically, the dual-spin satellite can be modeled 
as a system of two rigid bodies coupled together in such a manner that a 
rotation about an axis fixed in both bodies is the only relative degree 
of freedom. A general, exact analytical solution for the torque-free 
attitude motion of an arbitrary dual-spin satellite does not exist. 
However, Cochran, Shu and Rew, 13 building on previous work (Ref. 12, 
page 37) have obtained a complete, exact, analytical solution for the 
case in which: (1) one body is axisymmetric; (2) the axisymmetric body 
rotates relative to the other about the former’s axis of symmetry; (3) 
the other body is asymmetric; (4) the axis about which relative rotation 
occurs is the axis of major, or minor, moment of inertia of the asymmet- 
ric body; (5) the relative rotation is either free (no internal torque), 
or the relative spin rate is constant . A summary of the mathematical 
model of Ref . 13 is presented in Appendix B for completeness . 
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3.3 Multi-Rotor Satellite Model 

An exact analytical solution is also available for the torque-free 
attitude motion of a particular model of a satellite which contains 
spinning rotors; i.e., "momentum wheels," or "reaction wheels." In this 
solution, also due to Cochran and Shu, 14 the satellite physical model 
(see Fig. 4) is that of a gyrostat, consisting of two axisymmetric , 
constant-speed rotors, or momentum wheels, and a "carrier" rigid body 
which has a mass distribution such that the centroidal inertia ellipsoid 
of the system of rigid bodies is axisymmetric. 14 The model is general 
enough to represent any axisymmetric satellite that contains an 
arbitrary number of axisymmetric rotating components which together 
produce a resultant relative, or "internal," angular momentum vector 
that is not parallel to a principal axis of the system. As noted, the 
solution does, however, require constant internal angular momentum and 
an axisymmetric system . The solution is given in Appendix C. 




IV. SATELLITE MODEL FOR NUMERICAL SIMULATIONS 


4 .0 Introduction 

Many satellites contain variable-speed spinning rotors or other 
moving parts. Should an exact analytical solution for the torque-free 
attitude motion of such a satellite exist, it will undoubtedly be much 
more complicated than the single-rigid-body analytical solution 
discussed in Chapter 3 and given in Appendices A, B and C. Thus, for 
complex models, a practical alternative is to use numerical integration 
techniques to obtain even the torque-free, or "unperturbed" solutions. ^ 
Furthermore, to determine the motion of a satellite which is dynamically 
complex and/or is exposed to significant perturbing torques (for 
example, a LAM satellite in LEO), numerical integration of the equations 
of motion is the only approach generally applicable. In this Chapter, a 
physical model of a satellite is described. The corresponding 
mathematical model presented following the physical model has been 
incorporated in a digital simulation code along with models of 
gravity-gradient and aerodynamic torques. This code utilizes a 
fourth-order, Runge-Kutta, numerical integration algorithm to produce 
attitude motion time-histories. 

The purposes of the simulation code are: (1) to verify the 

correctness of the analytical solutions; (2) to determine the conditions 
under which torque-free solutions can be considered reasonably valid; 
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and (3) to produce attitude motion time-histories for cases in which 
analytical solutions are not known. 

4.1 Satellite Physical Model 

The "general" satellite model (see Fig. 5) consists of two rigid 
bodies, each of arbitrary mass distribution, which are coupled together 
so that the smaller body, or rotor, can be rotated with respect to the 
main body of the satellite, only about an axis which is fixed in the 
satellite and the rotor, and passes through the center of mass of the 
rotor. This physical model can be used to represent single-rigid-body 
and dual-spin satellites and also satellites which possess internal 
angular momentum due to constant-speed multiple, spinning components. 
Since the equations of motion for the general model are to be integrated 
numerically, environmental torques of some complexity may be included. 

Because the description of the attitude motions of satellites 
requires the use of reference frames and coordinate systems, the 
coordinate systems in Figs. 6, 7, and 8 are introduced. They are all 
dextral, orthogonal systems. The nonrotating EXYZ (inertial) 
coordinate system, depicted in Fig. 6, has its origin at the center of 
the earth and is considered inertially fixed. The unit vector triad 

AAA 

(I,J,K) is fixed to the EXYZ system. The coordinate system, Ex o y o Z Q , 
is called the "orbital plane system” and is defined by requiring that 
the XQ-axis lie along the line of nodes of the satellite's orbit, 
directed toward the ascending node, and that the z Q -axis be 
perpendicular to the orbital plane, collinear with the orbital areal 

A A A 

velocity vector. The unit vector triad is fixed to the 
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Fig. 7. Orbital Plane and Angular 
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Fig. 




. Angular Momentum and Body-Fixed Coordinate Systems . 
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Ex 0 y 0 z o coordinate system. An "orbiting coordinate system" C^nC which 
has its origin at the center of mass of the satellite is also shown in 
Fig. 6. The £-axis of this coordinate system is collinear with II, the 
position vector of C. The r)-axis lies in the orbital plane and is 
directed in the sense of increasing true anomaly, f. The unit vector 

AAA 

triad is attached to C The orientation of the orbiting 

frame is defined by using the standard angles: n, the longitude of the 

ascending node; i, the inclination; w, the argument of periapse; and f, 
the true anomaly. The notation, u”aj+f, is used for convenience. 

For cases in which the satellite's rotational angular momentum is 
large (HAM) compared to the magnitude of the external torque, the 
orientation of the satellite with respect to the inertial system is 
specified by first introducing the "angular momentum coordinate system," 
CxHyH z H> which has its zjj-axis collinear with H, the total angular 
momentum vector of the satellite due to rotation about its center of 
mass (see Fig. 7). The angle ’Pg is used to locate the xg-axis which 
lies along the intersection of the orbital plane and the plane which 
passes through C and is perpendicular to H. The angle between these 

A 

2 ,a^) is associated with the 
Cx HyH z H system. Finally, the orientation of the satellite-fixed, 
centroidal, principal system Cx^y^Zb with respect to Cxgyiizjj system is 
defined by using the Euler angles ij;, 0, and 4> (as shown in Fig. 8) in a 
3-1-3 rotation sequence. The angles ij>, 0, and <j> are identified as the 
classical angles of precession, nutation, and proper rotation, 

AAA 

respectively • 12 The unit vector triad is 

attached to Cx^y^z^. 


two planes is 0^ . A unit vector triad (a^ ,a 
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Alternatively, for the LAM cases in which the satellite's angular 
momentum vector is small in magnitude compared to that of the external 
torque, a "local vertical" system Cx v y v z v is used which has its z v -axis 
directed parallel to the negative £-axis of the C£n£ system and its 
x v -axls parallel to the orbital velocity vector. The orientation of the 
satellite body-fixed system Cx^y^^ with respect to the Cx v y v z v system 
is defined by the angles 63 , 82 , and 0^ (as shown in Fig. 9) in a 
3-2-1 rotation sequence. Angles 0 ^ , and 0^ correspond to roll, 
pitch, and yaw angles, respectively. 

4.2 Equations of Motion 

The equations governing the attitude dynamics of the "general" 
satellite model depicted in Fig. 5, may be derived using the angular 
momentum of the satellite, expressed in the satellite's body-fixed 
components; P a , the angular momentum of the rotor about the axis of 
rotation of the rotor with respect to the satellite; a, the angle of 
relative rotation; and four Euler parameters, gj, j=0,l,2,3. 

Let K denote the centroidal inertia matrix of the complete 
satellite, less the contribution by the rotor about the axis of rotation 
of the rotor with respect to the satellite. Then, the total angular 
momentum can be expressed as 

H = K u> + h , (4.2.1) 

where m is the absolute angular velocity of the satellite in the 
satellite's body-fixed components (b-basis) and where Ji , the internal 
angular momentum vector, is 
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h = J 


0 

0 

L V° 


(4.2.2) 


where J is the inertia matrix of the rotor about Cr, the center of mass 
of the rotor, in the b-basis; L^is transformation matrix from the 
r-basis (a coordinate system fixed in the rotor with its 3-axis along 
the axis of rotation of the rotor with respect to the satellite) to the 
b-basis; and D is the moment of inertia of the rotor about the rotor 
axis of rotation. 

The differential equation for the matrix, H, is 

H = H K -1 (H-h) + T , (4.2.3) 


where 


H 


0 

H 

2 

-H 


-H 


H 


z y 
0 -H 


H 


0 


(4.2.4) 


and ^ is the external torque expressed in the satellite’s body-fixed 
components. The momentum P a varies according to 


*C - (0 o 1)[H r/c (W +a)]+T (4.2.5) 

r ' R r r 

where is the angular momentum vector of the rotor about its own 

center of mass, expressed in the r-basis; u) and ^ are the angular 

r r 

velocities of the satellite and the rotor, respectively, also expressed 
in the r-basis; and T a is the torque on the rotor about the axis of 


relative rotation. 
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The Euler parameters 3 j , j=0,l,2,3, vary with time according to 

£^0 = " j | (4.2.6) 


and 




(4.2.7) 


where J = ( 3^ S^) 1 * (see, for example, Ref. 17, pp. 17-18 and 26). 
The direction cosine matrix, which defines the orientation of the 
Cx^y^z^ system with respect to the inertial coordinate system EXYZ may 
be expressed in terms of 3 ^ , j=l,2,3, as follows: 


*bl 


i - 2(3| + e§) 

2 <- hh - W 

2 <W w 


W 

1 - 2(02 + e|) 
2(S 2 8 3 - 6 t B 0 ) 


2< e i ®3- 8 2 8 o > 

2 <W B l®„> 

1 - 2(6 1 + B§) 


(4.2.8) 



V. TORQUES 


5.0 Introduction 

The environmental torques which act on artificial satellites are 
very small compared to, for example, the aerodynamic torques which act 
on aircraft. However, when even small torques act over significant 
periods of time, large changes in attitude motion can occur. As 
discussed in the chapter on ’’Categorization," if the time period of 
concern is relatively short, and if the magnitude of the satellite's 
rotational kinetic energy is relatively large compared to the work done 
on the satellite by its environment, the change iii a satellite's 
attitude will occur very slowly. A very important factor governing the 
effects of environmental torques on the attitude motion of satellites is 
the ratio of |t| , the magnitude of the perturbing torque to |Hj , the 
magnitude of the satellite's total angular momentum. When the ratio 
|T|/|H| is very small, say (5^(1 0” 1 * ) , as it would be for a HAM satellite, 
the short-term effects may be considered ’’minute." However, over long 
periods of time they can significantly perturb a HAM satellite's 
attitude motion. When |tJ/|h| is much larger, say^Cl), as it might be 
for a LAM satellite in LEO, the satellite's attitude motion will be very 
strongly perturbed, even in the short term. In some cases, tumbling 
motion which exhibits no well-defined pattern may occur. Under other 
circumstances a LAM satellite will tend to oscillate "about a position 
of relative equilibrium." 6 ’ 9 
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Generally, environmental torques fall into two main categories* 

The first category consists of environmental torques which are derivable 
from a potential function. The most studied member of this category is 
the gravity-gradient torque; however portions of the "aerodynamic” 
torque (due to interactions of the atmosphere and a satellite) may also 
be derivable from a potential function. These torques are therefore 
"conservative." However, when the motion of the center of mass of the 
satellite is specified, there is generally no integral of the motion 
corresponding to conservation of energy, but in the special case of a 
circular orbit, a Jacobi-type integral does exist if only "conservative” 
torques are present. The other category of torques is that of 
"dissipative," or "damping," torques. These nonconservative torques 
cause secular changes in the rotational kinetic energy of the satellite 
and its attitude motion. The more significant of these are: the 

dissipative parts of the aerodynamic torques; the torques due to the 
solar radiation pressure — direct from the sun, as well as reflected by 
the earth and its atmosphere; and torques resulting from the interaction 
of conducting parts of satellites with the earth’s magnetic field. The 
most dominant in orbits of altitude up to 300 km are the aerodynamic 
torques . 6 > 9 > 1 ® 

Many of the perturbation techniques employed in the treatment of 
satellite orbital motion can be employed analogously to study satellite 
attitude motion in the presence of perturbing environmental torques. 

Any further discussion, or reference, to environmental torques in this 
report will center on the two most dominant, gravity-gradient and 
aerodynamic . 
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5»1 Gravity-Gradient Torques 

According to Beletskii, 9 the primary effect of gravity-gradient 
torques "amounts to secular precession of the angular momentum vector 
around the normal to the orbital plane." He goes on to say that 
"periodic nutations of the angular momentum vector (with a period 
comparable to the satellite’s orbital period) are superimposed on this 
secular precession." 9 

Gravity-gradient torques are due to the fact that the earth's 
gravitational field is not uniform. This non-uniformity in both the 
magnitude and direction of the field over distances even as small as the 
dimensions of an orbiting artificial satellite result in a "gravity- 
gradient" potential over the satellite body. This gravity-gradient 
potential over the satellite results in a gravity-gradient torque about 
the satellite's center of mass. This gravity-gradient torque was first 
considered in the context of celestial mechanics in 1749 by d'Alembert 
and Euler, and later in 1780 by Lagrange . 1 1 • 1 8 

A number of factors are involved in determining the perturbation 
effect of gravity-gradient torque upon the satellite's attitude motion. 
If one assumes that the satellite is a single-rigid-body orbiting a 
spherical primary in a circular orbit, the problem is greatly 
simplified. With these assumptions, the gravity-gradient torque becomes 
a function of the distance from the center of mass of the satellite to 
the center of the earth; the values of the principal moments of inertia; 
and the orientation of the satellite's body-fixed axes with respect to 
the radius vector between the center of the earth and the satellite's 
center of mass. Using the above assumptions, the gravity-gradient 
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torques on an asymmetric satellite in a Newtonian gravitational field 
can be expressed in the form 9 6 1 9 


T 

"g 


3y E 


R 


(C B) c 23 c 33 
(A-C) c 33 c 13 


(B A) c 13 c 23 


(5.1.1) 


where R is the magnitude of R, the position vector from the center of 
the earth to the satellite's center of mass; A, B, and C are the 
principal centroidal moments of inertia; and the cj 3 , j=l,2,3, are the 
direction cosines of II in the satellite's body-fixed coordinate system. 
If the above simplifying assumptions are not made the expression for the 
gravity-gradient torque is much more complex. For example, if one 
assumes an "oblate earth," the expression for the x^-component of the 
gravity-gradient torque has the form 17 
3 U E 3p J r2 

\ = 13- <C-B) c 23 c 33 + -ssr* (C-B) [5(1-7 sin^c^ 

-10 sin X(c 23 c e33 + c 33 c e 23)“ 2c e23 c e33^ (5.1.2) 

where J 2 = 1.083xl0 -3 , \ is the latitude of the center of mass of the 
satellite, R e is the equatorial radius of the earth, and c e j 3 , j=l,2,3, 
are the direction cosines between the body-fixed axis and the earth's 
polar axis. Note that the additional gravity-gradient torque terms in 
Eq • (5.1.2) diminish at a rate of R -5 while the terms in Eq . (5.1.1) 
diminish at a rate of R -3 . For small satellites, the additional terms 
due to the asphericity of the primary can be neglected as they normally 
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result in less than a 1% difference, even for satellites in "low earth 
orbits." However, if one were to consider a larger satellite, such as 
the proposed Space Station, then the higher order terms would take on 
greater significance . 1 ? 

Another concern is that of more than one celestial body. In the 
earth, moon and sun system one must consider the gravity-gradient 
effects of the other celestial bodies in much the same manner that 
one considers their perturbing effects upon orbital trajectories. For 
near earth orbits, one may assume that the satellite is within the 
"sphere of influence" of the earth, and include only the perturbing 
effects of the earth. To illustrate this point, consider that for a 
geosynchronous orbit the gravity-gradient torque due to the moon is less 
than 0.0023% of that due to the earth, and the gravity-gradient torque 
at geosynchronous altitudes due to the sun is less than 0.0011% of that 
due to the earth . 1 f 

To summarize, for the purposes of the simulations used in this 
investigation, the following assumptions were found to give a 
sufficiently accurate representation of gravity-gradient torques on the 
satellites modeled: (1) the satellite mass distribution is that of a 
single, tri-inertial (or asymmetric), rigid body; (2) its center of mass 
is moving in a two-body orbit about a spherical primary; (3) the 
greatest dimension of the satellite is much smaller than the radius of 
the orbit of its center of mass; and (4) the orbital plane may be 
rotating in a prescribed manner. 
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5.2 Aerodynamic Torques 

As noted in the introduction to this chapter on environmental 
torques, there are a number of dissipative torques. The most dominant 
of these in the LEO altitude range are the aerodynamic torques. 

Although well outside the domain of "atmospheric vehicles," artificial 
earth satellites are by no means free of the effects of the atmosphere. 
This is demonstrated by the constant attention which those who track 
satellites pay combating the problem of orbital decay due to atmospheric 
drag. The treatment of aerodynamics at orbital altitudes as opposed to 
lower altitudes is different. At orbital altitudes the "mean free 
paths" between molecules are on the order of a kilometer or more, much 
larger than the dimensions of LEO satellites. This means that one need 
not consider collisions between molecules approaching the satellite’s 
surface with those leaving the surface. Thus, one can use a 
"free-molecular flow model" to approximate the "aerodynamics" of 
satellites. This treatment greatly simplifies the determination of the 
aerodynamic torques acting upon a satellite . 9 / 

There are four primary aerodynamic effects on satellites moving in 
the free-molecular flow regime. The first effect is similar to the 
"Weather Cock” effect which occurs in the lower portions of the 
atmosphere. This effect is due to the fact that the center of pressure 
and the center of mass of the satellite do not coincide. As a result of 
this, there is a restoring torque which tends to stabilize the satellite 
in the direction of the resultant of the orbital velocity of the 
satellite and the local velocity of the atmosphere due to the earth’s 
rotation. The second effect is due to the angular velocity of the 
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satellite about its own center of mass. The torques developed due to 
this effect are known as the "spin damping dissipative torques . ” 6 » 9 
Since even for high rates of satellite rotation the linear velocity due 
to rotation is much smaller than the satellited orbital velocity, these 
torques can be considered to vary linearly with angular velocity. The 
last two effects are small by comparison and are usually ignored. These 
are pressure gradient and the molecular thermal velocity effects. The 
pressure gradient effect is due to the density differential which exists 
over the satellite surface (being higher on the portion of the satellite 
facing the earth)* The contribution of molecular thermal velocity is 
generally neglected because it is much smaller than the orbital velocity 
of the satellite. 9 

For the purposes of this research effort, only the first two 
aerodynamic effects were considered to be significant. The model chosen 
to simulate the satellite in a free-molecular flow regime is a right- 
circular cylinder, a shape which is common among satellites. The 
coordinate axes chosen, as seen in Fig. 9, are very similar to those 
used by Etkin 20 to describe the motion of atmospheric vehicles. Note 
that the axis of symmetry has been designated the x^-axis in a 
^ X A^A Z A system. For the purposes of this work, it is assumed that 
the axes of C^x^yAZA are aligned with the axes of Cx^y^z^, but that C 
and Ca do not coincide. The CAX w y w z w , "relative wind" system, has its 
x w -axis along the orbital velocity vector and its z w -axis lies in the 
orbit plane and points (nominally) toward the center of the earth. 

The angles, a a and 0 a , between the relative wind system and the 
body-fixed system are the angle of attack and the sideslip angle, 
respectively, of the satellite (see Fig. 10). For the purpose of 
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integrating the force per unit area, the cylinder is divided into three 
surfaces, S 3 , S 2 and S 3 (see Fig. 11). Position vectors, r^ , r 2 , and 
r^ from the centroid of the cylinder to the surface elements dS^, dS 2 > 
and dS 3 , respectively, are used in deriving an expression for the torque 
about C^. The aerodynamic torque on the cylinder may be written as 
follows: 

T a - - 1/2 p c d [J r L x (n^VJV dS L + / r 2 x(n 2 -V)V dS 2 
S 1 S 2 

+ / r 3 x (n 3 .V)V dS 3 ] , (5.2.1) 

S 3 

where p is atmospheric density, is an accommodation constant, 

A A A 

n^ , n 2 , and n 3 are the normal unit vectors normal to each respective 
surface, V = + w x r ^ , j=l,2, or 3, where V c is the translational 

velocity of the satellite's center of mass,* and u> = (p q r)T. The 
results of the derivation of the aerodynamic torques on the right 
circular cylinder are given in Appendix D. 

It must be realized that any model of the aerodynamic torque on a 
satellite in orbital space is an approximation. Two factors which limit 
the accuracy of the modeling are variations in atmospheric density, and 
"shadowing,” or "blanketing,” of portions of the body by its other 
portions, for example, solar panels . 1 ^* 21 Since the purpose of the 
simulations contained in this thesis is to determine qualitatively the 


*Note that strictly speaking the velocity V_ should be used. 

However, this would only modify the results by changing the r^, 
j=l,2,3, slightly. 
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effects of environmental torques, a rather simple exponential 
atmospheric density model was adopted* Obviously, this model does not 
take into account such sources of variations in density as the diurnal 
cycles, orbital inclination (latitude variations at equal altitudes), 
and solar activity. Also, although a right-circular cylinder is a good 
approximation to the shape of many satellites, for satellites which have 
complex geometries (such as those with large sun seeking solar arrays) 
shadowing is a major problem. Such arrays can, at certain satellite 
attitudes, blanket large portions of the satellite surface resulting in 
torques greatly different from those on an unshadowed body . 21 However, 
the modeling of these complex satellite arrays for all possible aspects 
is extremely difficult and is beyond the scope of this effort. 



VI • SHORT-TERM APPLICABILITY OF TORQUE-FREE 
ANALYTICAL MODEL 

6.0 Introduction 

When a closed-form analytical solution exists for the torque-free 
motion of a particular satellite, it is generally the more efficient 
solution for short-term attitude motion simulations* Numerical 
solutions, although generally available, tend to require much more 
computer time than torque-free analytical solutions, especially in the 
case of satellites with high rates of rotation which necessitate the use 
of small integration step-sizes* However, when one attempts to 
incorporate the perturbing effects of environmental torques into the 
analytical approach, the complexity of the analytical solution increases 
significantly* In some cases, specifically, those in which the 
magnitude of the total angular momentum is much larger than the 
magnitude of the environmental torques acting upon the satellite, the 
torque-free analytical solutions are very accurate over short time 
periods* In other cases, for example, those in which the magnitude of 
the environmental torque is greater than the magnitude of the total 
angular momentum of the satellite, a numerical approach incorporating 
environmental torques is usually necessary* 

The purposes of this chapter are to investigate the fidelity and 
applicability of the analytical solutions discussed in Section 3.3 and 
to present some results of numerical simulations of LAM satellite atti- 
tude motion* First, to verify their fidelity, the analytical solutions 
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are compared with results from the numerical simulation code. Second, 
to determine the applicability of the torque-free analytical solutions, 
they are compared with numerical simulation results obtained using the 
same initial conditions, but including gravity-gradient torques. This 
comparison illustrates how, at lower angular momentum levels (as 
compared to the magnitude of the environmental torque), the fidelity of 
the torque-free analytical solutions degenerates rapidly, indicating the 
need to use either an analytical or a numerical solution which 
incorporates environmental torques. Third, to show that there are cases 
for which the analytical solutions cannot be expected to provide 
accurate results, some examples of simulated attitude time-histories of 
actual satellites in the LAM category are presented. 

6.1 Applicability of Analytical Solutions 

To verify the analytical solutions, the numerical simulation code 
was run with the data shown in Table 1. The data used in the analytical 
solutions is given in Table 2. Note that the moments of inertia in 
Table 1 differ slightly from those in Table 2. This is due to the fact 
that the numerical code uses a single movable rotor to simulate the two 
rotors of the analytical model. In each, the system is axisymmetric . 

As verification of the analytical solution, the results for the four 
cases are shown in Figs. 12, 13, 14, and 15. The variable "TAU" is a 
"no idimensional" time, defined as TAU = Ht/A. The numerical results are 
shown as triangles superimposed on the corresponding analytical 
time-histories (continuous curves) of the precession, nutation, and 
proper rotation angles, respectively. The first three cases (shown in 
Figures 12, 13, and 14) are for an oblate satellite with various initial 
angular velocities. The fourth case (Fig. 15) is for a rapidly spinning 


Table 1 Data for Numerical Verification of Torque-Free Analytical Simulation of Attitude 
Motion of Satellites with Asymmetric Internal Angular Momentum. 
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Orientation of Along y-axis In x-y plane Same as In x-y plane 

Rotor Spin Axis 71.56° from Case II 7.6° from 

x-axis x-axis 


Table 2 Data for Torque-Free Analytical Simulation of Attitude Motion 
of Satellites with Asymmetric Internal Angular Momentum. 
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1000*0 
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Fig. 12 a. Multi-Rotor Numerical Verification Case I, Precession 
Angle Time-History. 
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Fig. 12 b. Multi-Rotor Numerical Verification Case I, Nutation Angle 
Time-History. 
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Fig. 12 c. Multi-Rotor Numerical Verification Case I, Angle of Proper 
Rotation Time-History. 
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TAU 


Fig. 13 a. Multi-Rotor Numerical Verification Case II, Precession 
Angle Time-History. 
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Fig. 13 b 


Multi-Rotor Numerical Verification Case II, Nutation Angle 
Time-History. 
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TAU 


Fig. 14 a. Multi-Rotor Numerical Verification Case III, Precession 
Angle Time-History. 
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TAU 


Fig. 14 b. Multi-Rotor Numerical Verification Case III, Nutation Angle 
Time-History. 
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Fig. 14 c. Multi-Rotor Numerical Verification Case III, Angle of Proper 
Rotation Time-History. 
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TAU 


Fig. 15 a. Multi-Rotor Numerical Verification Case IV, Precession Angle 
Time-History. 
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prolate, dual-spin satellite with a misaligned rotor. As can be seen, 
the numerical and analytical results are, within the numerical 
precession of the calculations, the same. 

The numerical simulation code was modified to demonstrate the 
fidelity of the analytical solution for HAM cases in the presence of 
environmental torques. The modifications to the code included the 
incorporation of gravity-gradient torques, and a routine which would 
allow the sampling of the time-histories of the three attitude angles at 
a time many nutation periods in the future. The computer program was 
used to integrate the equations of motion over the intervening period, 
but output to data files for plotting was limited to the desired sample. 
The case which was run used the data from Case III in Table 1. The 
attitude motion of the satellite as perturbed by gravity-gradient torque 
was simulated for a time equal to a quarter of an orbit (1342.6 sec for 
an orbital altitude of 250 km). As can be seen in Fig. 16, the torque- 
free analytical time histories and the projected numerical solutions 
agree in amplitude very well. There is a “phase shift" which is prob- 
ably due to to the fact that the nutation period was not exactly equal 
to the number used to determine the beginning of the second sample time. 

To determine the magnitude of rotational angular momentum for which 
there would be significant disagreement between the two solutions, the 
initial angular momentum was successively reduced. To obtain the 
results shown in Fig. 17, the initial angular velocity and internal 
angular momentum components were each divided by 120. The results were 
chosen for illustrative purposes because they show a large difference in 
the amplitudes of the solutions, but the attitude motion structure has 
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Fig. 16 a. Fidelity of Multi-Rotor Results in the Presence of 
Gravity-Gradient Torque, HAM Case, Nutation Angle 
Time-History. 
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Fig. 16 b. Fidelity of Multi-Rotor Results in the Presence of 
Gravity-Gradient Torque, HAM Case, Angle of Proper 
Rotation Time History. 
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Fig. 17 a. Fidelity of Multi-Rotor Results in the Presence of 
Gravity-Gradient Torque, LAM Case, Nutation Angle 
Time-History. 
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Fig. 17 b. Fidelity of Multi-Rotor Results in the Presence of 
Gravity-Gradient Torque, LAM Case, Angle of Proper 
Rotation Time-History. 
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not yet broken down. That is, the three angles of precession, nutation, 
and proper rotation are still appropriate for describing the motion of 
the satellite. As was noted previously in Section 4.1, at low angular 
momentum levels, roll, pitch, and yaw angles are more appropriate for 
describing the satellite’s attitude motion because then 0, and <f> are 
not well defined. The amplitudes of the LAM (torqued) numerical results 
are significantly greater than the HAM torque-free analytical results. 
One could infer from these that the torque-free analytical solutions 
would be appropriate for short-term HAM satellite attitude motion 
simulations, and that the (torqued) numerical simulations more 
appropriate for short-term LAM satellite attitude motion simulations. 

6.2 Low Angular Momentum Examples 

The capability of the numerical simulation code when applied to low 
angular momentum (LAM) satellite simulations, was demonstrated by making 
eight simulation runs using data for three actual satellites. Four 
simulation runs were made using data for the Hubble Space Telescope 
(HST) . Then, using the data from the Earth Resources Satellite 
(LANDSAT) and the Advanced X-Ray Astrophysics Facility (AXAF), two 
computer runs were made for each of the two satellites. The satellite 
physical data for the HST satellite simulations can be found in Table 3. 
The corresponding orbital data can be found in Table 4. The data used 
for the LANDSAT and AXAF simulation runs found in Tables 5 and 6 is 
presented in the same manner. 

All of the results from the simulations are in the form of 
time-histories of selected parameters over five orbital periods. 

These parameters are: (1) the Euler parameters, 3 0 , , 62 » anc * 63; 


Table 3. Hubble Space Telescope (HST) Numerical Simulation, Satellite Physical Data. 



60 


Orientation of 

Rotor Spin Axis: Along xb-axis Along xb~axis Along Xb~axis Along xb“axis 



Table 4 Hubble Space Telescope (HST) Numerical Simulation Data. 
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Pitch 


Table 5 LANDSAT and AXAF Numerical Simulation, Satellite Physical Data. 
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Table 6 LANDSAT and AXAF Numerical Simulation Data 
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Pitch 
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(2) orientation of the satellite with respect to its local vertical 
coordinate system in terms of roll, pitch, and yaw angles; (3) the 
body-fixed components of the satellite's total angular momentum vector; 
(4) where applicable, the body-fixed components of the aerodynamic 
torque, and (5) again, where applicable, the components of the gravity- 
gradient torque, also with respect to the satellite's body-fixed axes. 

The four HST cases run include circular LEO, circular HEO, and 
elliptic LEO cases, all with environmental torques, and an elliptic LEO 
case without environmental torques. The results for the HST are shown 
in Figures 18, 19, 20, and 21. The LANDSAT cases are a circular LEO 
case and a circular HEO case. Both include only gravity-gradient 
torques. The two AXAF cases run are both circular LEO, one with 
environmental torques and one without. The results of the LANDSAT 
simulations are found in Figures 22 and 23. The results of the AXAF 
simulations are found in Figures 24 and 25. A more detailed analysis of 
the simulation results is made in the following paragraphs. 

6.2.1 HST Simulations 

In the first HST case, as with other HST cases, the satellite is 
initially "pitching down,” or rotating, about its negative y^-axis at 
orbital rate. The satellite also has some internal angular momentum due 
to a reaction wheel spinning about its xjj-axis. All of the HST cases 
begin with a 45° pitch angle with respect to the local horizontal. 

HST Case I 

The first figure for the HST (see Fig. 18 a) shows the 
time-histories of the Euler parameters, henceforth referred to as the 
"Beta plots." The individual Euler parameters oscillate in two modes; a 


65 


"low" frequency mode with a comparatively large amplitude and a "high** 
frequency mode with a much smaller amplitude. The time-histories of the 
roll, pitch, and yaw angles are presented in Fig. 18 b. Note that the 
pitch angle begins at 45° and oscillates at a relatively high frequency 
at progressively lower amplitude due to aerodynamic damping. After five 
orbits, the amplitude of the oscillation is approximately half its 
initial value. The roll angle shows relatively steady oscillations at 
low frequency, approximately equal to the orbital frequency and a 
smaller amplitude high-frequency mode. The yaw angle shows the same 
high frequency oscillation with fairly regular pulse-like changes in 
amplitude . 

The third figure (see Fig. 18 c) shows plots of the angular 
momentum components in which the damping of the component about the 
y^-axis is commensurate with the damping in pitch. The angular momentum 
component about the x^-axis remains nearly constant, and is due to the 
internal angular momentum and the lack of environmental torques about 
that axis. The x^-component is similar to the yaw oscillations in that 
it has relatively small amplitude pulses. 

The variations in the components of aerodynamic torque, shown in 
Fig. 18 d, are proportional to the corresponding angular momentum plots 
in all three components. In order to give the satellite model a margin 
of static stability, the center of gravity was placed 1/10 of the length 
of the satellite forward along the x^-axis from center of pressure (the 
centroid of the right-circular cylinder model) . This was also done for 


the AXAF satellite. 
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18 a. Time-Histories of Euler-Parameters for HST, 
Case I (250 km. Orbital Altitude). 
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Fig* 18 b. Time-Histories of Roll, Pitch, and Yaw Angles for HST, 
Case I (250 km* Orbital Altitude). 




Fig. 18 c. Time-Histories of Angular Momentum Components for HST, 
Case I (250 km* Orbital Altitude) * 
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Fig. 18 d. Time-Histories of Aerodynamic Torque Components for HST, 
Case I (250 km. Altitude). 
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Fig. 18 e. Time-Histories of Gravity-Gradient Torque Components for 
HST , Case I (250 km. Orbital Altitude). 
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The fifth, and last, set of plots for this case is given in 
Fig. 18 e. These are the time-histories of the components of the 
gravity-gradient torque. Here again the variations in the three 
components are proportional to those of the corresponding components of 
the vectors shown in the previous three figures. Notice, however, that 
the maximum amplitude of the oscillations in the gravity-gradient torque 
components is two orders of magnitude smaller than that for the 
aerodynamic torque . 

HST Case II 

Results for the second HST case were obtained using all of the same 
initial conditions as the first case except of orbital altitude. The 
first HST case was a circular orbit with an altitude of 250 km. In the 
second case, the orbit model was a circular orbit with an altitude of 
800 km. The Beta plots for this case (see Fig. 19 a) show somewhat 
similar low frequency/high amplitude oscillations of the parameter, but 
the high frequency/low amplitude oscillations have disappeared. The 
attitude angles (see Fig. 19 c) are very different from that for the 
first case. The pitch angle shows oscillations of similar amplitude 
initially, but does not exhibit the same aerodynamic damping. Also, the 
frequency of the oscillation in pitch is much lower for this case. The 
roll and yaw angles both oscillate at lower frequencies but at greater 
amplitudes than in the previous case. Also, note that roll angle 
oscillation is changing secularly in the negative direction. The 
angular momentum plot (see Fig. 19 c) shows that the x^-component 
remains constant. Note that the maximum amplitude the oscillations is 
an order of magnitude less than the first case. This is due to the 
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Fig. 19 a. Time-Histories of Euler-Parameters for HST, 
CASE II (800 km. Orbital Altitude). 





Fig. 19 b. Time-Histories of Roll, Pitch, and Yaw Angles for HST , 
Case II (800 km. Orbital Altitude). 


ANGULAR MOMENTUM COMPONENTS (KG-M** 2/SEC) *10* 

til 4 * A AA A ' A* 


74 



Fig. 19 c. Time-Histories of Angular Momentum Components for HST, 
Case II (800 km. Orbital Altitude) . 





Fig. 19 d. Time-Histories of Angular Momentum Components for HST, 
Case II (800 km. Orbital Altitude). 
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Fig. 19 e. Time-Histories of Gravity-Gradient Torque Components for 
HST, Case II (800 km. Orbital Altitude). 
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lower orbital rate and low level of excitation by aerodynamic torque. 

The y^-component is varying secularly in the positive sense while 
oscillating slowly with amplitude which decreases and then increases. 

The z^-component exhibits what appear to be long-period variations. 

While changing in this manner, the z^-component is also oscillating much 
like the y^-component . The aerodynamic torque components (see Fig. 19 
d) are several orders of magnitude smaller than they were in the first 
case; in fact, they are smaller than the gravity-gradient torque 
components (see Fig. 19 e) . This is due to the fact that the model of 
the density of the atmosphere decreases exponentially with altitude, 
whereas the gravity-gradient torque is proportional to the inverse cube 
of orbital radius. 

HST Case III 

The third HST case is a low, slightly elliptic orbit (e = .01). As 
a result of the ellipticity, the orbital altitude ranges from a perigee 
altitude of 184 km to an apogee altitude of 316 km. The Beta plots in 
Fig. 20 a exhibit two frequency modes seen in Case I. But the 
amplitudes of the lower frequency/high amplitude modes show some 
parameters with growing amplitudes and others with diminishing 
amplitudes. The attitude plot for this case (see Fig. 20 b) shows that 
the roll angle oscillates in two frequency modes as in the first case, 
but it is also changing secularly in a negative sense. The pitch angle 
time history contains a high frequency oscillation (as in the first 
case) and is decreasing in amplitude as time progresses, but there are 
amplitude pulses which initially coincide with apogee passage. These 
pulses most likely result from the lower level of pitch damping due to 



Fig. 20 a. Time-Histories of Euler-Parameters for HST, 

Case III (184 km. x 316 km. Orbital Altitude). 
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Fig. 20 b. Time-Histories of Roll, Pitch, and Yaw Angles for HST , 
Case III (184 km. x 316 km. Orbital Altitude). 
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Fig. 20 d. Time-Histories of Aerodynamic Torque Components for HST , 
Case III (184 km. x 316 km. Orbital Altitude). 
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Fig. 20 e. Time-Histories of Gravity-Gradient Torque Components for 
HST, Case III (184 km. x 316 km. Orbital Altitude). 
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lower atmospheric density at the apogee altitude. The yaw angle 
oscillations have approximately the same frequency as the oscillation in 
pitch. However, the yaw angle oscillations begin with zero amplitude 
and then increase. Initially, the yaw amplitude pulses coincide with 
the pitch amplitude pulses in frequency. As time progresses and the 
amplitude of the yaw angle pulses grows larger than the diminishing 
pitch angle pulses, it also becomes out of phase. The yj,- and 
z ^-components of angular momentum in Fig. 20 c behave very similarly to 
the pitch and yaw angles, respectively, of the attitude plot. The 
x^-component of angular momentum remains nearly constant due to the 
internal angular momentum of the satellite and the fact that there is no 
aerodynamic torque (Fig. 20 e) about the x^-axis . The y^-component of 
aerodynamic torque decreases in amplitude with periodic pulses which 
coincide with perigee passage. The z^-component of aerodynamic torque 
remains at a relatively low amplitude with pulses also coincident with 
the perigee y^-component pulses. The amplitude of the y^-and 
z^-components of the gravity-gradient torque pulse not so much as a 
result of the changing orbital altitude as with the changing of 
orientation due to the rolling motion of the satellite. That is, when 
the y^-axis is pointed more towards the center of the earth, the 
zjj-component of gravity-gradient torque increases in magnitude. The 
converse is also true. 

HST Case IV 

This case has the same initial conditions as Case III. However, in 
this simulation the environmental torque subroutines were "disconnected" 
so that the motion is torque-free. The Beta plots (Fig. 21 a) show only 
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Fig. 21 a. Time-Histories of Euler-Parameters for HST, 

Case IV (184 km. x 316 km. Orbital Altitude without 
Environmental Torques). 
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Fig. 21 b. Time-Histories of Roll, Pitch, and Yaw Angles for HST, 
Case IV (184 km. x 316 km. Orbital Altitude without 
Environmental Torques). 
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Fig. 21 c. Time-Histories of Angular Momentum Components for HST, 
Case IV (184 km. x 316 km. Orbital Altitude without 
Environmental Torques) . 
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the expected low frequency/high amplitude oscillations. The attitude 
angLes (Fig. 2L b) show the same secular variation in roll angle and 
almost the same for low frequency oscillation as in Case III, but the 
high frequency oscillation has disappeared. The pitch angle has a low 
frequency/low amplitude oscillation and slight secular departure in the 
negative direction. The yaw angle shows a slightly higher amplitude and 
a frequency similar to that of pitch angle. The yaw angle does not 
appear to change secularly. The angular momentum (Fig. 21 c) has a 
constant x^-component , and the y^-and z^-components exhibit single- 
frequency, low-amplitude oscillations, but are out of phase with each 
other. The low amplitude is due to the low angular rates which are 
unexcited by environmental torques. 

6.2.2 LANDSAT Simulations 

The LANDSAT attitude motion simulations differ from the HST simula- 
tions in a few respects. First, the satellite physical data indicates 
that the satellite’s moments of inertia are an order of magnitude 
smaller than the HST. Second, the LANDSAT model possesses no internal 
angular momentum as did the HST. Third, the initial orientation of 
LANDSAT model with respect to the local vertical differs from that of 
the HST. The HST simulations began with an initial pitch angle of 45°, 
but the LANDSAT simulations begin with a 10° pitch angle and a -80° yaw 
angle. Fourth, the initial angular rate of the HST was a negative pitch 
about the y^-axis at orbital rate, while the initial angular rate of the 
LANDSAT is ncoslO 0 about the x^-axis and nsinlO 0 about the z^-axis , 
where n is the orbital mean motion for a given orbital radius. Fifth, 


88 


the orbits which are used for the LANDSAT simulations are both inclined 
at 98.2° as compared to the 28.5° inclination used for the HST cases. 

The first LANDSAT case is for a circular orbit with an orbital 
altitude of 352 km. The second LANDSAT case is also for a circular 
orbit, but the orbital altitude is 705 km. In both of the LANDSAT 
cases, only gravity-gradient torques were simulated. 

Case I 

In the 352 km orbital altitude case, the Beta plots (Fig. 22 a) are 
irregular oscillations which show strong orbital coupling. These 
oscillations are most similar to the HST Case IV Beta plot (Fig. 21 a), 
but are somewhat more irregular. The attitude angle plots (Fig. 22 b) 
shows that roll angle increases secularly with some irregular 
"oscillations." The pitch angle time-history also shows similar 
irregular "oscillations," but no secular departure. The yaw angle 
time-history contains higher amplitude irregularities and is also not 
changing secularly. The angular momentum plot (Fig. 22 c) indicates 
that this case is a much lower angular momentum case than the HST cases. 
Note that the scale used in the plots is two orders of magnitude smaller 
than that for the previous HST angular momentum component plots. The 
xjj-component remains nearly constant and the y^~and z^-components have 
low frequency/low amplitude oscillations and no apparent secular change. 
The gravity-gradient torque component plot (Fig. 22 c) shows that the 
gravity-gradient torque is an order of magnitude less than it was for 
the HST cases. This is due to the fact that the LANDSAT satellite is 


much smaller than the HST . 
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Fig. 22 a. Time-Histories of Euler-Parameters for LANDSAT, 

Case I (352 km. Orbital Altitude, Gravity-Gradient Torque 
Only) . 
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Fig. 22 b. Time-Histories of Roll, Pitch, and Yaw Angles for LANDSAT, 
Case I (352 km. Orbital Altitude, Gravity-Gradient 
Torque Only) . 
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Fig. 22 c. Time-Histories of Angular Momentum Components for LANDSAT, 
Case I (352 km. Orbital Altitude, Gravity-Gradient 
Torque Only) . 
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Fig. 22 d. 



Time-Histories of Gravity-Gradient Torque Components for 
LANDSAT, Case I (352 km. Orbital Altitude, Gravity- 
Gradient Torque Only) . 




Fig. 23 a. Time-Histories of Euler Parameters for LANDSAT, 

Case II (705 km. Orbital Altitude, Gravity-Gradient 
Torque Only) . 
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Fig. 23 b. Time-Histories of Roll, Pitch, and Yaw Angles for LANDSAT, 
Case II (705 km. Orbital altitude, Gravity-Gradient 
Torque Only). 
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Fig. 23 c. 


Time— Histories of Angular Momentum Components for LANDSAT , 
Case II (705 km. Orbital Altitude, Gravity-Gradient 
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Fig. 23 d. Time— Histories of Gravity— Gradient Torque Components for 
LANDSAT , Case II (705 km. Orbital Altitude, Gravity- 
Gradient Torque Only) . 



97 


Case II 

The results of the second case are nearly identical to that of the 
first case. The Beta plot (Fig. 23 a) are almost identical to that for 
the first case. The attitude angles for both cases (Figs. 22 b and 23 
b) are also almost identical. The angular momentum component plots 
(Figs. 22 c and 23 c) are slightly different. The orbital mean motion 
in this case is very slightly lower due to the fact that the orbital 
altitude is 705 km. As a result of the slightly lower angular velocity 
magnitude, the angular momentum components have slightly lower 
amplitudes. Again in the gravity-gradient torque plot (Fig. 23 d) , the 
amplitudes are slightly smaller due to the higher orbital altitude. 

6.2.3 AXAF Simulations 

This is the final satellite considered in the numerical 
simulations. The AXAF satellite is more similar to the HST than the 
LANDSAT. Geometrically, the sizes of the AXAF and the HST are nearly 
identical. However, the moments of inertia of the AXAF are slightly 
more than twice as large as those of the HST. The AXAF model includes 
no internal angular momentum. The initial orier Nation of the satellite 
model has a 45° pitch angle with respect to the local horizontal as does 
the HST. Also, the AXAF model has a 45° yaw in the negative sense. The 
initial angular rates of the AXAF satellite model are l.OxlCF 3 rad/sec 
about the x^-axis and l.OxlO” 4 rad/sec about the y^-axis . The orbital 
inclination is the same as for the HST simulations, 28.5°. Both of the 
AXAF simulations assume a circular orbit with an altitude of 400 km. In 
the first case, both of the environmental torques were simulated. In 
the second case, no environmental torques were included. 






TIME (SEC) *10 


Fig. 24 b. Time-Histories of Roll, Pitch, and Yaw Angles for AXAF , 
Case I (400 km. Orbital Altitude) . 


I 
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Fig. 24 c. Time-Histories of Angular Momentum Components for AXAF, 
Case I (400 km. Orbital Altitude). 
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Fig. 24 d 


Time-Histories of Aerodynamic Torque Components for AXAF , 
Case I (400 km. Orbital Altitude). 
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Fig. 24 e. Time-Histories of Gravity-Gradient Torque Components for 
AXAF, Case I (400 km. Orbital Altitude). 
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Case I 

The Beta plots (Fig. 24 a) show the reappearance of the two 
frequency mode oscillations seen in the first HST (LEO) case. Here, the 
high frequency/low amplitude oscillations have a lower frequency and are 

more irregular than those in the first HST case. This is due to the 

high frequency/ low amplitude oscillations have a lower frequency and are 

more irregular than those in the first HST case. This is due to the 

intermediate altitude of this case (400 km) , which is between the first 
HST case (250 km) and the second HST case (800 km) where the high 
frequency oscillations have disappeared. The roll angle in Fig. 24 b 
behaves much differently than in previous cases. The roll angle 
increases to approximately 650°, then decreases to around 200°. The 
pitch and yaw angles appear to be almost periodic with no secular 
variations. They behave much like the pitch angle in the second HST 
case (see Fig. 19 b) . The components of angular momentum (Fig. 24 c) 
also behave in a manner which suggests that this case represents a 
'* transition" between the first and second HST cases. Note that the 
maximum amplitude of the angular momentum is an order of magnitude 
greater in this case. There is no apparent damping effect due to 
aerodynamic torque on any of the components. Note that the x^-component 
is not constant as it has been in previous cases. The reason for this 
can be seen in the gravity-gradient torque component plot (Fig. 24 e) . 
The x^-component of gravity-gradient torque is not constant because 
there is a large difference between the moments of inertia Iyy and I zz . 
The aerodynamic torque components (Fig. 24 d) are essentially 
proportional to the angular momentum components. Note that the maximum 
amplitude of aerodynamic torque lies somewhere between those for the 
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first and second HST cases. Also, note that the maximum amplitude of 
the gravity-gradient torques Is approximately twice that of the first 
HST case. This is due, as we stated previously, to the fact that there 
are greater differences in the values of the three principal moments of 
inertia. 

Case II 

This final numerical simulation case has the AXAF model with 

\.-s 

initial conditions identical to the first AXAF case. The differences; 
between the two are a result of the absence of environmental torques in 
the second case. 

The Beta plots (Fig. 25a) show only the single low frequency/high 
amplitude oscillations. The attitude angle plots (Fig. 25 b) show that 
all three angles "oscillate” at a relatively low frequencies with 
significant irregularities. The final plot of this section, the angular 
momentum plots (Fig. 25 c) , are similar to the corresponding plots for 
the second HST case. That is, the x^-component remains constant, while 
the other two exhibit smooth oscillations of the same frequency and 
amplitude and are slightly out of phase. 

6.2.4 Summary of Low Angular Momentum Examples 

The results of the eight cases given in this section indicate that 
the time-histories of some of the variables used to describe the 
attitude motion do not exhibit well-defined patterns. In some cases, 
aerodynamic and gravity-gradient torques perturb the attitude motion of 
low angular momentum satellites to such an extent that existing 
torque-free analytical solutions are not valid for any significant 
period of time. 




Fig. 25 a. 


Time-Histories of Euler-Parameters 
Case II (400 km. Orbital Altitude, 
Torques) . 


for AXAF , 

No Environmental 



I 



ATTITUDE ANGLES (DEG) 
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Fig. 25 b. Time-Histories of Roll, Pitch, and Yaw Angles for AXAF, 
Case II (400 km. Orbital Altitude, No Environmental 
Torques) . 


ANGULAR MOMENTUM COMPONENTS (KG-M**2/SEC) *10’ 





Fig. 25 c. Time-Histories of Angular Momentum Components for AXAF , 
Case II (400 km. Orbital Altitude, No Environmental 
Torques) . 



VII. PREDICTION OF LONG-TERM MOTION 


7 .0 Rationale 

As was stated previously, there will likely be a time lapse of 
several months between a satellite's "failure" and a recovery attempt. 
Therefore, a need exists for an accurate long-term attitude motion 
simulation for mission feasibility studies and planning. This chapter 
focuses on this need. As with the short-term attitude motion 
simulations, there are two approaches to the problem of perturbed 
long-term attitude motion. These are the General Perturbation Methods 
(GPM) , and the Special Perturbation Methods (SPM) . To apply the usual 
GPM, one needs an analytical solution to the unperturbed problem. 
Examples of unperturbed solutions include two-body orbital motion and 
"free-Eulerian" rotational motion. Tractable analytical expressions for 
perturbations are also required. The SPM most commonly applied are that 
of numerical methods. Examples of numerical perturbation methods 
applied to orbital motion include the work of Cowell and Encke. 24 
Applications to attitude motion include the work of Kraige and 
Junkins. 25 The short-term motion simulation, described in Chapter IV, 
also uses numerical perturbation methods. 

Both GPM and SPM have shortcomings which hamper their application 
to attitude dynamics problems. The application of the GPM is difficult 
or perhaps impossible if the satellite is dynamically complex. On 
the other hand, the numerical approach, although generally applicable. 
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often requires large amounts of computer time, especially when higher 
angular rates necessitate the use of very small integration step sizes. 
For these reasons, a "hybrid" method was developed and labeled the 
"Semi-Analytic method." This method incorporates the use of analytical 
averaging concepts from GPM, and numerical integration techniques from 
SPM to obtain "averaged equations of motion" which govern the long-term 
attitude motion of certain classes of satellites. 

t 

7.1 Averaging Method 

In order to apply the General Method of Averaging to the equations 
of motion to obtain the "averaged equations of motion," new variables 
must be introduced. These variables are elements of x and y, where 

x = (H Hcos0 H ¥ h ), (7.1.1) 

and 

I m (4> * M) , (7.1.2) 

where M is the mean anomaly. The elements of x are the "slow" variables 
and the elements of y. are the "fast" variables. The time derivatives of 
£ and _y take the forms , 

X = e fj^x* y) (7.1.3) 

and 

i ■ z) + y) , (7.1.4) 

where g Q (x, y) is the "unperturbed" time derivative of y, e f^(x, y) 
and e y) are the "perturbations" due to external torques on the 

satellite, and e is the usual small parameter introduced for 
perturbation analyses . 
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Here, it will be assumed that the functions of f^(x, y) and 
gj( x , y) > j = 0,1, are periodic with period 2tt in the elements of y« 
Since the elements of y_ are Euler angles , this will be the case for a 
variety of torques . 

A new set of variables (x, y) is defined by 


x = x + e u 1 (x, y) + e 2 y 2 (x, & + *** 

and 


(7.1.5) 


Y = y + e (x, y) + e 2 v 2 (x, y) + (7.1.6) 

where the functions Uj and v^, j = 1,2, are to be chosen to make the 
differential equations for x and y simpler than those for x and y* 
Equations (7.1.5) and (7.1.6) may be differentiated with respect to time 
to get 


x = x + e 


Hi + 


*2 + 


(7.1.7) 


and 

i = Y + e Yi + e2 £2 + ••• (7.1.8) 

The respective right-hand sides of Equations (7.1.3) and (7.1.7), and 
(7.1.4) with (7.1.8) may be equated to provide 


x + e u x + e 2 v 2 + ... = e f-^x, Z) 


(7.1.9) 


and 


Y + e v x + e 2 v 2 + ... = g Q (x, Z) + e Z> 


(7.1.10) 


where the arguments of and have been dropped for brevity. 

The aim of the averaging method is to obtain averaged equations of 
motion which do not contain the "fast variables." It is assumed that 
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the Uj and Vj are & (1) (or smaller), and the coefficients of equal 
powers of e in Eqs . (7.1.5) and (7.1.6) are grouped, then the 
averaged equations of motion may be obtained in the forms , 


x = e f 1 (x, 2 ) + e ^(x. X) 


(7.1.11) 


and 


i = s 0 ($> i) + e S!<i» i) 


(7.1.12) 


where an overbar over a function indicates its "average" value. This 
requires that the "periodic perturbations" be chosen such that 

H 1 = fj^x, 2) - f x (X» jr) (7.1.13) 

and 

ii = Si(i> i) - Si(i» £) + Os 0 /3x) 0 ^ 

+ (3S 0 /3Z) 0 Hi » (7.1.14) 

where ( ) 0 indicates evaluation of the function ( ) at e = 0. The 

function can be found by putting the zeroth-order solutions for 

x and y into and integrating with respect to time. Then, may be 

found by solving Equation (7.1.14), which is a linear differential 
equation if the zeroth-order solutions for x and y are used. Here, 

"average” is defined using the zeroth-order solutions for x and y. For 

a f^*ee, single-body satellite, the angles ^ and are ordinarily 

monotonic, while terms like sin0flsin<j) have zero averages. Other terms, 

for example, cos0, usually have non-zero torque-free average values. 
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In torque-free motion, when the satellite contains spinning rotors ({> may 
be bounded. A further explanation of this behavior is given in Ref. 15. 

7 .2 Equations of Motion 

The equations of motion for this semi-analytic method were derived 
from expressions for the angles ©h, <j>, tJj, 0 and |h| , the magnitude 
of the total angular momentum vector. The angles were described in 
Section 4. The angles and ©h describe the long-term motion of the 
angular momentum coordinate system with respect to the orbital plane 
coordinate system (see Fig. 7). The angles ij;, 0, and <j> describe the 
short-term motion of the satellite's body-fixed coordinate system with 
respect to the angular momentum coordinate system (see Fig. 8). In 
vector notation, the total rotational angular momentum vector can be 
written as 

H = I • m + h , (7.2.1) 

where I is the principal, centroidal inertia dyadic, matrix, w is the 
angular velocity vector, and li is the internal angular momentum vector. 
In matrix notation, this expression may be written as 

(H x H y V T = I a + (h x h y h z ) T (7.2.2) 

An expression for the angular velocity vector can be derived from 
Eq. (7.2.2). It is found that 

ui = I -1 (H - h) . (7.2.3) 

The time rate of change of the angular momentum vector, H, written in 
vector form is 

H = H + a> x H 


9 


(7.2.4) 
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where 

H = (ft ft H ) (a. a_ a_) T 
x y z -1 -2 -3 

In matrix notation, the external torque expression is 
T = ft - H to , 


(7.2.5) 


( 7 . 2 . 6 ) 


where 


H = 


0 

H 

z 

-H 


-H 


H 


z y 

0 -H 

x 

H 0 


(7.2.7) 


and 


ft = (ft ft ft ) T . 

x y z 

H can be written in the satellite's body-fixed system, or in the angular 
momentum system. In the satellite system, it is written 

ft = H to + T . (7.2.8) 

Alternatively, by substituting in Equation (7.2.3), one finds 

ft = H I -1 (H - h) + T . (7.2.9) 


In the angular momentum system, 


= 5h % + • 


(7.2.10) 


The 


where ^=(00 H) T , H = |h| and T R = (T T y T z ) T . 

a H H 

angular velocity of the angular momentum system, u^, includes the 
secular effects of the regression of the line of nodes of the 
satellite's orbit. Explicitly, 
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4 ! 


ft si st 


H 


ft [si c0 }1 ct H + ci sO H ] 
ft[ci c0 H - si s0 H c'Fj,] 



h 1 


+ 

*H Sln0 H 

= 


• 



» H co S 0 H - 



H H 


0) 


4! 


03 


r H 


0) 


(7.2.11) 


The relationship between o) and is co = o^ 4- X , where 


c<j> 

s9s<{) 

0 

-S(J> 

S0C(j) 

0 

0 

c9 

1 


1 

•CD 

1 


• 

¥ 


1 

• -©- 


(7.2.12) 


is the angular velocity of the satellite with respect to the angular 
momentum system. 

The equations of motion for the short-term variables include 
expressions for a, b, and ij. The parameters a and b are used to 
define 0, the nutation angle, where a and b are defined as 


and 


a = H sinQ 


b = H cos0 . 


(7.2 13a) 
(7.2.13b) 


The components of the angular momentum vector can then be found from the 
relations 



H 


H 


= a sin<j), 
= a cos<J), 


(7.2.14a) 

(7.2.14b) 
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and 

H z = b. (7.2.14c) 

The time derivatives of a and b, which may be integrated to find the 
time-history of the nutation angle, are 

a = b{[a sin<j) - h ) cos<j >]/ A - [(a cos<{> - h ) sin(j>]/B} 
x y 

+ sin<j> + Ty cos (7.2.15) 

b = [(B-A)/AB] a sin<f> cos<f> 4* (lWA)a cos<J> 

- (h y /B) a sin<j> + T z> (7.2.16) 

where A, B, and C are the principal centroidal moments of inertia, 
and where T z is the zy-component of external torque. From Eqs . 

(7.2.14), (7.2.15) and (7.2.16), one may obtain the following equations 
for <J> and if>: 

l - (b-h z )/C-b[B sin 2 <ji + A cos 2 <(,]/AB 
+ (b/a)(A h cos<j> + B h sin<j>)/AB 

y x 

+ (T /a) cosij> - (T /a)sin<j> (7.2.17) 

x y 

and 

J - H[B sin 2 <f> 4- A cos 2 <f)J/AB - {[Ah cos <J> 4- Bh sin<{>]/AB 

y x 

+ [(T /H) cos\jj + (T /H) sin\{>] cos9}/sin9 

*H y H 

- (T /H) cot0 - Q sini cosf u /sln0 . (7.2.18) 

X u n rl ri 

H. 

The time derivatives of the long-term variables fn, 0H> and H are 

= - Q sin i cosf„ cotG.. - Q cos i + (T /H)sinQ, I , (7.2.19) 

rl H X_ T rl 


= - n sin i sin^u - T /H 
H y H 


(7.2.20) 
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and 

ft = T , (7.2.21) 

H 

where T , T , and T are the components of external torque in the 
X H y H Z H 

angular momentum system. 

7 >3 Results 

To demonstrate the capability of the semi-analytic simulation 
method, two examples are given. The first example 15 involves an 
axisymmetric satellite which contains two axisymmetric rotors on the 
satellite's x^-and y^-axes , respectively. The second example concerns 
the Combined Release and Radiation Effects Satellite (CRRES). 26 In this 
section, a brief description of the computer code designed to use the 
semi-analytic method is first given. Next, the results of the two 
examples are presented and finally, a summary of these results is made. 

For descriptive purposes the semi-analytic computer code is 
broken down into four parts: the initialization process; the short-term 

stage, or Stage 1; the intermediate stage, or Stage 2; and the long-term 
stage, or Stage 3. 

In the initialization process, I, the principal centroidal inertia 
matrix; o) Q , the initial angular velocity vector; h, the internal 
angular momentum vector; and the initial values for the angles 
and 0 , are entered. From I, u) , and h, the total angular momentum 
vector, H, are then used to find the initial values of the angles 
0 and <j). The initial gravity-gradient torque is also calculated to find 
the initial values for Stage 1. 
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The six equations of motion for the Stage 1 integration are those 
for: <J>; \p; 4^; P^, the momentum conjugate to \J;, where = H; P , the 

momentum conjugate to <|>, where P = b; and P , the momentum conjugate 

♦ H . . . , . 

to , where P w = H cos0 H . The expressions for <j> , \J>, H, b, and 


H 


H 


H 1 


0 can be found in the previous section. The nutation period is 
H 

determined either by using an analytical solution, as in the first 
example, or numerically, as in the second example. The six equations of 
motion are then integrated over one nutation period to find average 
values of P , P t , P, and \b, which are to be used as initial conditions 
for the ’’averaged equations of motion” for Stage 2. This completes the 
short-term stage, or Stage 1. 

The intermediate stage, or Stage 2, uses these "averaged equations 
of motion." These four averaged equations of motion are integrated 
numerically, over an orbital period. The integration step size is much 
larger (two orders of magnitude) than the step size used in Stage 1. 

The results are doubly averaged values of P^ and ^ which are used then 

averaged again over the orbital period to obtain as the initial values 

• • 

= = 

for the "doubley averaged equations of motion" for P and V . 

In the final stage, or Stage 3, the two doubly averaged equations 
of motion are integrated over a specified "long" period of time using 
integration step sizes on the order of one-half day. 

Exrmple 1 

Data for the first example is given in Table 7. Using the 


analytical solution, a nutation period of 4.6 seconds was obtained and 
the mean time rates of change of $ and were found to be approximately 
1 .367 rad/sec and 2.132 rad/sec, respectively. Plots of <j>, 0, and are 
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shown in Figure 26 for the Stage 1 portion of the semi-analytical 

solution process. The mean value of 9 during the short-term motion is 

approximately 6.267 deg. For Stage 2, the numerically determined motion 

in ijl, P , and is shown in Figure 27. Results for Stage 3, or 

long-term motion is shown in Fig. 28 in the form of a 0„ vs. 4' plot. 

ri n 

The analytical solution to the averaged equations yields the same 
information. 


Table 7 . Data for Semi -Analytic Example 1 


A = B = 400 kg-m 2 
C = 200 kg-m 2 
h x = 20 kg-m 2 /s 

hy = 0 

h z = 150 kg-m 2 /s 
0 ^( 0 ) = 0.1 rad/s 
ujy(0) = 0.001 rad/s 
io z (0) = 3.5 rad/s 


i|>(0) = 0 
'fjj(O) = 0 

0 H (0) = 80 deg 

R => 6778.27 km 
e = 0 

i = 28.5 deg 
ft * -6 deg/day 
n = 3.998 deg/min 


As a partial check on the accuracy of the semi-analytical theory, 
the full equations of motion were integrated numerically using a step 
size of 0.2 seconds. Results using a fourth-order Runge-Kutta algorithm 
on a Harris 800 minicomputer over 29,164 seconds of CPU time are 
compared with the semi-analytical results in Fig. 29. The small initial 
difference in ¥ and is due to using the average value from Stage 2 
of the semi-analytical approach. The CPU time for the application of 
the semi-analytic theory is around 12 seconds. Clearly, good accuracy 
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Fig. 29 


Comparison of Direct Integration Results and Semi 
Analytic Results for Semi-Analytic Example 1. 
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can be maintained with tremendous savings of computer time using the 
semi-analytic, theory. 

As a second example, long-term attitude motion of the CRRES 
satellite was simulated using the serai— analytic method. The data for 
this simulation is presented in Table 8. Using a numerical solution, a 
nutation period of 3.9 seconds was determined. Plots of <)>, 0, and \J* are 
shown in Fig. 30 for the Stage 1 portion of the semi-analytic process. 
The mean value of 0 during the short-term motion is approximately 4.58 
degrees. For Stage 2, the numerically determined motion in ij>, P , and 
and is shown in Fig. 31. Results from Stage 3, for long-term motion, 

n 

are shown in Fig . 32 in the form of a Gjj vs f ^ plot . 


Table 8. Data for Semi-Analytic Example 2 (CRRES) 


A = 
B = 
C = 

h x = 
h v = 


‘z 

o)x(0) 

U>y(0) 


h, = 


: ( 0 ) = 


2263.13 kg-m 2 
1917.5 kg-m 2 
3719.65 kg-m 2 

0. 

0. 

0. 

0.15 rad/s 

0 . 

1.0472 rad/s 


f(0) = 0 

*h(°> = o 

0h(O) = 5 deg 

R = 7378.27 km 
e = 0 

i = 28.5 deg 
Q = —6 deg/day 
n = 3.4246 deg/min 


NUTATION ANGLE (DEG) 


TIME (SEC) 

(•) 


(SEC) 

<b) 



TIME (SEC] 


Semi -Analytic Example 2 
Time-Histories, Stage 1 
(b) Angle of Proper Rot. 


irt-Term Motion, 
on Angle; 
ecession Angle. 
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The CPU time required for the second example was 9.22 seconds. The 
integration time step size was 0.1 sec for Stage 1, 20 sec for Stage 2 
and 1.2 day for Stage 3. 

A semi-analytic method for predicting the long-term attitude motion 
of uncontrolled satellites has been described. The method can be used 
effectively when the square of the rotational angular momentum of the 
satellite under consideration is much larger than the maximum magnitude 
of the torque multiplied by the largest moment of inertia of the 
satellite. The method is applicable in both cases in which the 
unperturbed solution is known and in cases for which it must be 
determined numerically. The problem of determining the attitude motion 
of a satellite with constant internal angular momentum and subjected to 
gravity-gradient torque was used in the first example. The problem of 
determining the attitude motion of a satellite which has no internal 
angular momentum, but is spinning quite rapidly and is subjected to 
gravity-gradient torque was considered in the second example. The speed 
with which solutions can be obtained using the method, relative to that 
of straightforward numerical integration is of the order of 2430 to 1. 
Such time savings are of considerable significance. 


VIII. CONCLUSIONS AND RECOMMENDATIONS 


8.1 Conclusions 

Much of the emphasis of this investigation has been placed on 
analytical solutions for the attitude motion of uncontrolled satellites . 
Analytical solutions to the torque-free equations of rotational motion 
for three satellite models were verified using numerical integration. 

The majority of this verification was done using a digital simulation 
code which contains a rather general dynamic model of a two-body 
satellite. The applicability of one of the torque-free analytical 
solutions for the attitude motion of satellites with asymmetric internal 
angular momentum was tested using a model for the gravity-gradient 
torque. It was determined that for high angular momentum cases the 
fidelity of the torque-free analytical solution was very good. However, 
for cases in which the satellite is modeled as slowly tumbling in the 
presence of environmental torques, the torque-free analytical solution 
does not accurately model the motion. Numerical results were obtained 
using physical characteristics of actual satellites for cases of low 
angular momentum in the presence of environmental torques . 

The investigation of long-term attitude motion emphasized the 
semi-analytic method. Good results with a substantial savings in CPU 
time were obtained using this method. It can be utilized in cases where 
the unperturbed analytical solution is known, and also those in which 
the unperturbed solution must be determined numerically. 
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8.2 Recommendations 

It is recommended that for short-term modeling of satellite 
attitude motion, for use in capture and retrieval simulators, the 
“driver" for dynamically simple, high angular satellites should be a 
“torque-free," closed-form analytical solution if available. This 
recommendation is based upon the “computational-speed advantage" of 
closed-form, analytical solutions. If, on the other hand, the satellite 
in question has such complex dynamics that an analytical solution is not 
available and/or is a low angular momentum satellite which is moving in 
the presence of significant environmental torques, a perturbed 
analytical solution or a numerical solution is required. Numerical 
integration of the attitude equations becomes more efficient in the low 
angular momentum regime because larger integration step sizes may be 
employed at lower angular rates. 

No analytical solutions for short-term prediction of the attitude 
motion of low angular momentum satellites are included in this report. 

It is recommended that such solutions be developed. These solutions 
would be very useful in satellite capture and retrieval operations. 

The recommendation for long-term prediction of satellite attitude 
motion is to use a semi-analytic method. Such a method should provide 
results which are sufficiently accurate for the purposes of predicting 
the state of the satellite's attitude motion sufficiently far in advance 
for capture and retrieval mission planning. 

The semi -analytic method presented herein does not address the 
problem of internal energy dissipation. Therefore, it is recommended 
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that "energy-sink" models be incorporated in the present serai-analytic 
method* 

Finally, for low earth orbit cases, due account must be taken of 
the aerodynamic torque. To do this precisely would require that the 
aerodynamics of each target in low earth orbit be modeled* It is 
therefore recommended that methods be developed to produce quickly 
attitude time histories of fairly general satellites perturbed by 


aerodynamic torques 
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A.l Introduction 

Many satellites may be considered to be "single bodies” which are 
fairly rigid. The simplest physical model for such a satellite which 
takes into account the possibility of mass asymmetries is an asymmetric, 
or ”tri-inertial rigid body. Although the solution for the torque- 
free rotational motion of such a body is well known 1 1 because it forms 
the basis for perturbed motion analyses and, principally, because it 
does provide a good approximation to the motion of many satellites over 
short time periods, a summary of the solution is included in this 
report. This summary closely follows that given in Ref. 27. 


A. 2 Mathematical Model 

Let the principal, centroidal moments of inertia of the body be A, 
B and C where either A>B>C or C>B>A. Let the motion of principal axes 
x b> yb and z b> measured with respect to a (fixed) rotational angular 
momentum coordinate system which has its zjj-axis collinear with Jfl, the 
angular momentum of the body about its center of mass C. Let the 
principal axes components of oj , the angular velocity of the body, be 
o>y , a^, and <o z . Also, let \J; , 9 and $ denote Euler angles such that 

= 


0 ) 


y 


\p cose sin<f> + 0 cos(f) 


(A. 2.1) 


a> z = cos9 cos<J) - 9 sin<() 
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Note that the Euler angle sequence used here is a 3-2-1 sequence . 

The equations of rotational motion for the body are 

Am + (C-B)oj a) = 0 
x y z 

Bm + (A-C) a) u = 0 (A. 2. 2) 

y z x 

Co) + (B-A)oj to = 0 
z x x y 

Also, the angular momentum can be expressed as either 

AAA 

H = Au> b 1 + Boo b„ + Cu) b_ , (A. 2. 3a) 

“ x “i y “Z z *"j 

or 

A A A 

H = H(-sin9 b^ + cos9 sin<(> + cos9 cos<}i b^), (A. 2. 3b) 

A 

where b., j=l,2,3, are unit vectors attached to the Cx.y.z. system. 

“ J d d b 

Equations (A. 2.1) have two immediate integrals, 

Ao) 2 + Bo) 2 + Co ) 2 = 2 °f (A. 2. 4b) 

x y z x 

and 

A 2 u| + B 2 a£ + C 2 u£ = H 2 , (A. 2. 4b) 

where ZF is the rotational kinetic energy. These integrals may be used 
in Eqs. (A. 2.1) to obtain 2 ^ the equation, 

l 2 = X 2 (l - k 2 sin 2 ?), (A. 2. 5) 

where the new variable £ = sin -1 (oo y /Q) ,Q = (d/e) 1/2 , d = 2 C tT -H 2 , 
e = B(C-B) , X 2 = [ (C-B) (H 2 - 2 A^)]/(ABC) and 
k 2 = [ (B-A) (2C Z) — H 2 ) ] / [ (C-B) (H 2 - 2 A^)] . 

It follows from Eq . (A. 2. 5) that 


sin£ = sn u 


(A. 2. 6) 
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where sn u is one of the Jacobian elliptic functions, u = Xt - v, and 
v is the value of u at t=0. Hence, 

ay = Q sn u . (A. 2. 7) 

For C>B>A, the integrals (A. 2. 4) amd (A. 2. 5) and Eq. (A. 2. 7) 
may be used to get 


uy = P cn u 


(A.2.8) 


and 


0 ) = R dn u , 

z 


(A. 2. 9) 

where P 2 = (2C fT - H 2 )/ [A(C-A) ] , R 2 = (H 2 - 2A?T)/ [C(C-A) ] and cn u and 
dn u are Jacobian elliptic functions. 

By using Eqs . (A. 2. 3) and defining p = AP/H, q = BQ/H and 
r = CR/H, the following relationships may be obtained from Eqs. (A. 2. 7) 
through (A. 2. 9): 

-sinG = p cn u , 

cosG sin<() - q sn u (A. 2. 10) 

and 

co s 0 cos<j> * r dn u . 

The angle of precession can be obtained from 


$ = h{ 1 + [ (C-B)/B] sin 2 <f>}/C , 


(A. 2. 11) 


which is derivable from Eqs. (A. 2.1) and (A. 2. 2). By using 
tan<b = (B/C) (q/r)(sn u/dn u) in Eq . (A. 2. 11), it is found that 


ip = i|> o + / (H/C)[l + g 2 sn 2 u/(l + a 2 sn 2 u)]dt , (A. 2. 12) 


where 


g 2 = a 2 (C-A)/A 
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and 


o 2 = [A(C-B) ] / [C(B-A) ] k 2 . 

The Integral in (A. 2. 12) is a form of elliptic integrals of the 


third kind and may be evaluated using methods described in Ref . 28 
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B.l Introduction 

A "dual-spin" satellite is usually a satellite which is composed of 
a nonrotating, or "despun,” "platform" and a rapidly rotating "rotor." 
More generally, a dual-spin satellite is one which has two parts which 
spin at different rates. The latter description is used for the 
purposes of this appendix. Additional requirements on the model of a 
dual-spin satellite adopted here are explained in the following# 

A solution to the equations of torque-free attitude motion of a 
particular dual-spin satellite physical model is given in this appendix. 
It is taken from Ref. 13. The physical model consists of a rigid, 
axisymmetric body and a rigid asymmetric body (see Fig. 3). Either can 
represent the platform (rotor). The axisymmetric body rotates with 
respect to the asymmetric body about its axis of symmetry which is 
parallel to one of the axes of extremal moment of inertia of the asym- 
metric body. For the solution given, the torque between the rotor and 
platform is assumed to be zero. However, the form of the solution is 
the same if the relative spin rate of the rotor is constant. 

B.2 Equations of Motion 

Let A*, B and C represent the principal, centroidal moments of 
inertia of the satellite about the x^-, y^- and z^-axes, respectively, 
which are fixed in the asymmetric body. Also, let denote the moment 
of inertia of the axisymmetric body about its symmetry axis, which is 
assumed to be parallel to the x^-axis. Furthermore, let A = A* - B| . 

The components of the rotational angular momentum are 

H = Aco + P (B.2. la) 

x x a 7 


B (0y 


(B.2. lb) 
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and 

H z = Cu> z , (B .2.1c) 

where u^, ay and u> z are the components of angular velocity of the 
asymmetric body and 

P a = B l (a) x + (B .2.2) 

is the angular momentum of the axisymmetric body about its symmetry 
axis . 

In Eq. (B.2.2), Si is the x-axis component of the angular velocity 
of the axisymmetric body with respect to the other body. 

Because the motion is torque-free, the angular momentum system 
CxRyH z H ( see Fig* B.l) is fixed and the Euler angles ip, 9 and <|> are 
defined so that 

H x = H cos0 , (B.2.3a) 

H y = H sin0 sin<t> (B.2.3b) 

and 

H z = H sin9 cos<{> , (B.2.3c) 

where H = (H x 2 + Hy 2 + H z 2 ) i/2 . Note that in this appendix and Appendix 
C, the XR-axis is aligned with H and the Euler angle sequence is 1-2-1. 

The equations of motion and the kinematic equations for \J> and a, 
the angle of relative rotation of the axisymmetric body, are 

\ = [(B-C)/(BC)]H y H z , (B .2 .4a) 

Hy = — {[(A-C)/(AC)]H x + P a /A}H z , (B.2.4b) 

fi z = {[(A-B)/(AB)]H x + P a /A}H y , 


(B .2 .4c) 
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II 

o 

(B.2.4d) 

a 

* = H(H y 2 /B + H z 2 /C)(H 2 -H x 2 ) 

(B.2.4e) 

a = P (A+B.)/(AB ) - H /A . 
a i ix 

(B.2.4f) 

Solutions for ip, 0, <p and a 

Equations (B.2.4) have the ’'immediate" integrals 13 

H 2 = H 2 +H 2 + H 2 , constant , 
X y z 

(B .3.1) 

P = constant 
a 

(B.3.2) 

2T = (H -P ) 2 /A + H 2 /B + H 2 /C, constant, 
x a y z 

(B.3.3) 


These three integrals may be used to write 

f 2 (H x ) = H 2 /{(b(A-C)]/[A(B-C)]} - [H x + C P^A-C)]* 

+ H 2 [A(B-C) ] / [B(A-C) ] - [H + C P /(A-C)) 2 (B.3.4a) 

y o x a 

and 


f 3 (H x ) = H z 2 /{[C(A-B)]/[A(B-C)]} = [H x + B P a /(A-B)] 2 

+ H 2 [A(B-C) [ / [C(A-B) ] - [H + B P /(A-B)] 2 , (B.3.4b) 
Z o o a 

where a subscript o denotes an initial value. 

Hence, Hy and H z are functions of H x , and H x 2 may be written in the 

form, 


fi x 2 = [(A-C)(A-B)/(ABC)]f 2 (H x )f 3 (H x ) . 


(B .3.5) 
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The solution to Eq. (B.3.5) depends upon whether the roots of the 
equations f2 “ 0 and f'3 = 0 are real or complex. For physical Ly 
realizable motion, the roots of f 2 = 0 "lust be real. The roots of 
f 3 = 0 may, however, be either real or complex. Hence, two cases are 
possible . 

Case 1 . Four Real Roots , a>b>c>d 
In this case , 

H x = (D^ + D 2 sn 2 u)/(D.j + sn 2 u) , (B.3.6) 

where sn u is a Jacobian elliptic function of modulus 

k = {[(a-b)(c-d)]/[(a-c)(b-d )]} l/2 (B.3.7) 

and argument , 

u = At - u 0 , (B.3.8) 

in which 

X = {[(A-B)(A-C)/(BC)] l/ 2 [(a-c)(b-d)] l/ 2 }/2A . (B.3.9) 

The forms of the D-:, j=l,2,3,4, depend upon H x . It turns out that 

J o 

H x is never between b and c for real motion. For b H x _< a, 
o o 

= a(b-d), T >2 - d(a-b) , D 3 = b-d and D 4 = a-b. For d _< H x _< c, 

o 

= d(a-c), D 2 = a(c-d), D 3 - a-c and D 4 = c-d. 

Case 2. Two Real Roots 

If only the roots of f 2 = 0 are real, the solution for H x has the 

form, 

H x = (C 1 + C 2 cn u)/(C 3 + cn u) . 


(B .3 .10) 
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Here, cn u is a Jacobian elliptic function of modulus 


k = { [ (a-b ) 2 - (a* - B*)2]/(4 a*g*)} 1/2 . 


(B.3.11) 


where a and b are the real roots (a>b) and a* and g* are the real and 
imaginary parts, respectively, of the complex root c. In the argument, 
u = Xt + u 0 , 


X = [(A-C)(A-B)/(A2BC)] 1 / 2 (a*B *) l/2 


(B.3.12) 


The constants C j , j-1,2,3,4, are defined as ■ a 6 * + ba*, 

C 2 = ba* - ag*, C 3 =■ a* + g* and C 4 = a* - g*. 

Solutions to the equations for if>, <}>, and a may be obtained for each 
of the two cases above. For Case 1, the solution for ip is 


where 


2 u 

<P “ l t> . / du/(l-a, 2 sn 2 u) + \Jj , 
J" 0 ' u o ' 


b o ^ C 1 + C 2 C 4 + C 3 C 6^ X 


“ C3C7 / X 


b 2 = c 2 c 5 /X 


(B .3 .13) 


(B .3.14) 


c t = H/A > 


c 2 = [2AT-(H+P a )2]/(2A) 
c 3 = [2AT-(H-P a )2]/(2A) • c 4>6 = D 4 /(D 4 H ± D 2 ) 

C 5,7 = (D 3 D 2 " D 1 D 4 )/[(D 3 H ± ± D 2 )] 


a Q = ° 


c ^ 2 = - (d 4 h - d 2 )/(d 3 h - d 3 ) 

c ^ 2 = - (D 4 H + D 2 )/(D 3 H + D x ) 


(B .3 .15) 


(B .3.16) 
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The integrals which appear in the solution for tp are elliptic integrals 
of the third kind (see Ref. 28, pp. 232-237), 

For Case 2, the solution for ip is 
2 


u 


4) = £ d. / du/(l - g. 2 cn u) + tp Q , 

i=0 J u - 

J A 


(B .3.17) 


where the dj and gj have the same forms as bj and otj but replaces 
in the definitions. Also, the modulus and parameter X have the required 
forms. The integrals in Eq. (B.3.17) are also elliptic integrals of the 
third kind (see Ref. 28, pp. 215 and 232-237). 

The solution for <j> has the form, 


<p = £ P i / du/(l - a. 2 sn 2 u) + d Q , 

i=0 J u J 

J o 

for case 1. The constants pj are defined as 
Po = [ p „/ A " C 3 C 6 + C ? C J/ X » 


(B .3.18) 


- C^Cy/X 


(B.3.19) 


and 


p 2 ~ C 2 C 5^ X 


For Case 2, the solution for $ is analogous to the solution for ip 
in that case. 

Solutions for a are similar to those for ip and <p, but are simpler 
because only one elliptic integral is involved. For Case 1, 

a = a Q (u - u q ) + / du/(l - a^ 2 sn 2 u) , 


(B .3 .20) 
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where 

a o = t p a ( A+B l )/B l - D 2 /D 4 ]/(AA), 

a l = (D 2 D 3 “ D i D 4> /(A D 3 D 4 *) (B .3.21) 

and 

a 3^ = ~ ° 4^3 * 

The solution for a for Case 2 can be obtained in the manner described 
above for and <(>. 
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C.l Introduction 

The model presented in Appendix B cannot be used to represent a 
satellite which contains a part (or parts) which is rotating about an 
axis which is not parallel to a principal axis of the carrier body. A 
model which can be used to represent such a satellite, if the system is 
axisymmetric , is presented in this appendix. In fact, the satellite 
model can be used to represent a satellite which contains several 
constant-speed, rigid "rotors.” Hence, the label, "Multi-Rotor 
Satellite." 

The physical model is depicted in Fig. 4. It consists of a 
torque-free carrier body and one, or more, constant-speed, axisymmetric, 
rigid rotors. These rotors are arranged so that their resultant angular 
momentum due to rotation relative to the carrier body (internal angular 
momentum) is not parallel to any axis, except possibly the y^-axis of 
the system Cx^y^Zb fixed in the carrier body. The x^-axis of the system 
is an axis of symmetry. Hence, the y^-axis can always be chosen such 
that the z^-component of the internal angular momentum is zero. The 
angular momentum system CxflyH z H is fixed because there are no external 
torques . 

Equations of Motion 

The mathematical model presented here is based, for the most part, 
on Ref. 14. Let H x , H y and H z denote the x^-, y^-and z^-components , 
respectively, of the total angular momentum about C. Also, let h x and 
hy denote the x^-and y^-coraponents , respectively, of the internal 
angular momentum and let A and C be the principal centroidal moments of 
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inertia of the system of carrier body and rotors, with A about the 
x^-axis . Then , the equations of rotational motion may be put into the 
forms , 14 

fl « H h /C , (C.2.1a) 

fi y = -{[(A-C)/(AC)]H x + h x /A}H z (C.2.1b) 

and 

ft z = {[(A-C)/(AC)]H x + h x /A}H y - H x h y /C . (C.2.1c) 

The corresponding kinematic equations are 

♦ - (H x - h x )/A - (H x /C)[l - h y H y /(H2-H x 2)] (C.2.2a) 

and 


4> = (H/C) [1 - h H / (H 2 -H 2 ) ] . (C .2 .2b) 

y y x 

Here, <|> and are the angle of proper rotation and the precession angle 
of the body, respectively (see Fig. B.l in previous appendix). 

The nutation angle 9 is given by 

0 = cos" 1 (H x /H) , (C.2.3) 

where H = (H x 2 + H y 2 + H z 2 ) 1/2 . 

The following expressions for the angular velocity components 
o) x , u)y y and (d z may also be required at times: 

u> x = (H x - h x )/A (C .2 .4a) 


Wy = (H y - h y )/c 

" H z /C 

Furthermore, the angles 9, <j> and can be used to write 

H = H cos 0 
x 

H = H sin© sind. 

y 


(C .2.4b) 
(C.2.4c) 

(C.2.5a) 

(C.2.5b) 
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and 


H z = H sin9 cos(|> , (C.2.5c) 

Equations (C.2.1) admit from "first" and "second" integrals 
H = constant and 


(H -h ) 2 /A + (H -h ) 2 /C + H 2 /C = 2T, constant 

xx y y 6 

Equation (C.2.6) may be rewritten (by using H»H = constant) 

H = aH 2 + 2bH + c , 
y x x 1 


(C.2.6) 


in the form, 
(C.2.7) 


where 


a = (C-A)/(2Ah y ) , 
b = - h x C/(2Ah y ) 

and 

c = [(C/A)h x 2 + h y 2 + H 2 - 2CT] / (2h y ) . 

Equation (C.2.7) represents a family of parabolic cylinders 
momentum space. The terminus of the vector II describes, on 
H 2 = H x 2 + Hy 2 + H z 2 , a curve which is an intersection of a 
cylinder with the sphere. 

These integrals may be used in Eq . (C.2.1a) to get the 


(C .2.8a) 
(C .2.8b) 


(C .2 .8c) 

in angular 
the sphere 
parabolic 


equation 


y ■ <y c > 2 y< H *> > 


where 


H z 2 ( H x> = a 2 [-H x ^ + c 


- H 3 
3 x 


+ c 2 H x 2 + c x H x c Q ] 


(C.2.9) 


(C .2 .10) 


c 3 = - 4 h x C/(A-C) (C .2.1la) 

c 2 = - 4(A 2 h y 2 + C 2 h x 2 )/(A-C) 2 

+ 2 [Ch 2 + Ah 2 + A(H 2 - 2CT)]/(A-C) 
x y 


(C.2.11b) 
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c L = 4[C 2 h x 3 + ACh x h y 2 + ACh x (H 2 - 2CT)]/(A-C) 2 (C.2.11c) 


c = A 2 {4h 2 H 2 - [ (C/A)h 2 + h 2 + H 2 - 2CT J 2 } / (A-C) 2 . 
o y x y 

(C.2.11d) 

C .3 Solutions for 6, <j> and i|< 

The general solution to Eq. (C.2.9) can be expressed in terras of 
elliptic functions of the time. The form of the solution depends on the 
number of real roots of the quartic equation, H 2 (H ) = 0. 

Z X 


Case 1. Four Real Roots 

In this case, if the roots are H Xj , 
k>j, then 


j=l ,2,3,4, and H x < H x for 
j k 


H x = (D x + D 2 sn 2 u)/(D 3 + D 4 sn 2 u) , (C.3.1) 

where sn u is Jacobi's sine amplitude function and the Dj are determined 

by H x , the initial value of H x . 
o 


If H > H > H , then 
X 1 - X o X 2 


D. = H (H - H ); 
1 x i x 2 x 4 


D„ = H (H - H ); 

2 V X 1 X 2 


If H x > H x > H x , then 
2 o 3 


D = H (H - H ); 
1 x 2 X 1 x 3 


D„ = H (H - H ); 
2 X 1 x 2 X 3 



(C.3.2) 


(C.3.3) 
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Finally, if H > H > H , then 
X 3 X o X 4 


D, = H (H - H ); D, = H - H 
1 x^ 3 x^ Xg 


D 0 = H (H - H ); D. = H - H 
2 X 1 x 3 x 4 4 x 3 x 4 


(C.3.4) 


The argument of the elliptic functions is 


u = At 4- u 


o y 


where 


X 2 = [ah /(2C)][(H - H )(H - H )] 

y x 2 x 4 x : x 3 


(C.3.5) 

(C.3.6) 


k 2 = {[(H - H ) (H - H )]/[(H - H )(H - H )]} . 

X 1 x 2 3 x 4 A 1 3 2 4 

(C.3.7) 

By using the solution for H x in Eq • (C»2.2a), and considerable 
algebraic manipulation, one may put the equation for <J> into the form, 

3 u 

= £ P. / du/(l - a. sn 2 u) , (C.3.8) 

j=0 J u J 

J A 


where a Q = 0, 

a x = -(D 4 /D 3 )d 9 /d 7 

02 = -(D 4 /D 3 )dg/d 6 
a 3 2 = " D 4^ D 3 

P q = [ (d 12 + d 4 dg)d 1 - d 5 d 9 d 2 + d 13 D 2 /D 4 ]/X (C.3.9) 

P 1 = d 5 d ll //X 
P 2 = d 4 d 10 /X 

P 3 = ~ D 2 D 3 )/(D 3 D 4 )]/X . 
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In writing Eqs. (C.3.7), use was made of the defined constants, 
d x = D 4 (D 3 H + D 1 )/[D 3 (D 4 H + d 2 )] 

d 2 = d 4 (d 3 h - D 1 )/[D 3 (D 4 H - d 2 )] 

d 3 = H(1 + ah y )/C 

d 4 = -h y (aH 2 - 2bH + c)/(2C) 

d 5 = -h y (aH 2 + 2bH + c)/(2C) 

d 6 = V< D 4 H + D 2> 

d 7 = D 6 /(D 4 H -D 2 ) (C .3.10) 

d 8 = V< D 3 H + D l> 

d 9 = D 3 /(D 3 H - D x ) 

d 10 = (°3 D 2 - D 4 D 1 )/[(D 4 H + D 2 )(D 3 H + D 1 )] 
d ll = (D 3 D 2 " D 4 D 1 )/[(D 4 H " D 2 )(D 3 H " D l> ] 

d 12 = "(V A + 2bh y /C > 

d 13 = 1/A - (1 + ah y )/C . 


The integrals in Eq . (C.3.18) are elliptic integrals of the third 
kind. A set of subroutines has been developed to evaluate these. 

By using a similar procedure , it may be shown that 


♦ - I 


jlo Qj l 1 " a i snZu + 




(C .3 .11) 


where 


Q o - (d 3 + d 4 d 6 + d 5 d ? )/X 

Qi = - p i 


(C .3 .12) 
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and 


^2 = 


Case 2. Two Real Roots 

When the quartic equation has two real roots, H x ^ > H x ^ » and two 
complex roots, y ± Si , the solution for H x has the form, 

H x = (C 1 + C 2 cn u ) / ( C 3 + c 4 cn u ), (C .3.13) 


where, if. 


and 


then 


r = [ (H - y) 2 + 6 2 ] L ' /2 

X 1 


s = [(H - y) 2 + 5 2 ] 1/2 , 

X 2 


C, = H s + H r , 
1 X 1 X 2 


C 0 = H r - H s , 
2 x x^ 


= r + s 


and 


C 4 = r - s . 


(C .3 .14a) 
(C .3 .14b) 


(C .3 .15a) 
(C .3 .15b) 
(C .3 .15c) 

(C.3.15d) 


The modulus of cn u is given by 

k = {[(H - H . ) 2 - (r - s) 2 ]/(4rs)}l/2 

X 1 x 2 

As in Case 1, u has the form, 
u = At + u Q . 

But, for this case, 

X = ah y [rs] 1/2 /C. 


(C .3.16) 


(C .3 .17) 


(C .3 .18) 
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Regarding the forms of the solutions for <J> and the only way in 
which this case differs from Case 1 is that the Cj are used in place of 
the Dj and the elliptic function cn u replaces sn 2 u. The elliptic 
integrals are still of the third kind (see Ref. 28, pp. 215, 232-237). 


APPENDIX D 


GRAVITY-GRADIENT AND AERODYNAMIC TORQUE MATHEMATICAL MODELS 
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D . 1 Gravity-Gradient Torque Mathematical Models 

Two mathematical models were derived during this investigation to 
model the gravity-gradient torque experienced by a satellite orbiting 
the earth. The first model was used for the short-term numerical 
simulation program. The second was used for the long-term semi-analytic 
simulation program. 

Short-Term Mathematical Model 

1 ' 1 ! | 

The unit vector e^ describes the attitude relationship between the 
satellite's body-fixed axes and the z v -axis of the local vertical 
system. Recall that the z v -axis always points toward the center of the 
earth. Thus, the third column of the L^ v transformation matrix provides 
the unit vector of interest in the matrix form, 



(D.1.1) 


This expression for e may then be used to complete the expression. 


T 

“8 



e 

R3 " 8 


I 


e 

-g 


(D.1.2) 


where is the gravitational constant of the earth, R is the distance 
from the center of the earth to the center of mass of the satellite and 
I is the matrix of the satellite. 

Long-Term Mathematical Model 

The model of the perturbing effects of gravity-gradient torque on 
the satellite attitude motion used in the semi-analytic simulation is 
based on a rigid axisymmetric body with centroidal moments of inertia 
A=B*C. The components of the gravity-gradient torque can be derived by 


160 


taking partial derivatives of the expression for gravity-gradient 
potential , 


V g = [3p E /(2R 3 )](C-A) cos 2 x , (D.1.3) 

where the angle x is the angle between the position vector of the 
satellite's center of mass and the body-fixed Zj,-axis (the axis of 
symmetry). In terms of the angles 0, On, i|>, and U • fn “ u » which are 

defined in the body of the report, the cos 2 x may be expressed as 

follows: 

cos 2 x = ^- [1 + c 2 0g + c 2 0 - 3c 2 0 H c 2 0] 

+ -^ [1 - c 2 0jj - 3 c 2 0 h c 2 0] cos 2 y 
+ -y 1 + c 2 ©^ "b c 2 0 — c 2 0] cos 2tp 

- [(1 + 2c0 h +c 2 0 h ) s 2 0] cos 2(i[»-p) 

+ y [(1 + 2 c 0 h ) s 2 0] cos 2(t|rfp) 

+ [s0gC0 H s0 c0] cos i|» 

+ j [(1 - c0 H )s0 H c0 s 0 ] cos ((|»-2p) 

- j 1(1 + c0 H )s0 H c 0 s0 ] cos (i|H-2p). (D.1.4) 

D.2 Aerodynamic Torque Mathematical Model 

The results of the derivation of the aerodynamic torques on the 
right-circular cylinder are as follows. The total aerodynamic torque 
acting on the right circular cylinder is a sum of the torques 
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T , T , and T contributed by the surfaces S. , S 0 , and S 

-“ci 2 "~a 3 i ^ 

respectively (as shown in Fig. 11). That is 


T = T + T + T 
-a - &1 -a 2 -a 3 


(D.2.1) 


The torque contributed by the front endcap, S]_, is 


T = (T T T ) T 

_a l a l a l a l 

x y z 

where 

\ ■}» c d r c“ ! 2UP - " r - V< »- 

X 


(D.2.2) 


T 3l = " i p c d r c 4 U[(r c 2 + 2x i 2)q ' 2Wx l ] * 

y 

and 

T a = - i P c d r c 4 U t (r c 2 + 2x l 2) r + 2VX 1 ] (D.2.3) 

where r c is the radius of the right circular cylinder, U, V, W are the 

x b“> Yb~> and z^-coraponents of the "relative wind,” where 

U = |V | cosa cosg , V = |v | sine , and W = |V | sina cose > and p, q, 
c a a ca caa 

and r are the components of angular velocity about x^-, y^-, and 
z^-axes, respectively, and x^ is the distance from the centroid of the 
cylinder to center of the front end cap along the x^-axis. The 
aeroiynamic torque contributed by the sides of the cylinder, surface 
S 2 , is 



(D.2.4) 
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where 

T 

a 2 

x 

T 

a 2 

y 


" I p c d r c { v [^vZ+W*y"i '(d^-dj 2 ) - fu r c L] 

- p[f» r (d,2-d 0 2)] - q [2 VW_ ( 2 L + d 3+^3)] 

4 C 1 2 3 /TvW) c 1 2 

+ r[|(/Tv^+wZT + — )(r 2 l + d ,3 + d„ 3 ) 

3 /'(V 2 +W 2 ) c 1 2 

“ £ u r c (d ^2- d ^2 ) ) ] } . (D.2.5) 

Here, L is the length of the cylinder, dj is the distance along the 
x^-axis between the front endcap and the center of gravity, and d 2 is 
the distance from the rear endcap to the center of gravity. Finally, 
the aerodynamic torque contributed by the rear endcap, S 3 , is 


and 


= - P c d r c 3 L /(V^+W*) , 


- - \ p c d r c {W(f U r. L - /TvWT (dj2 - d 2 2)] 


- p[|-V r c (d 1 2-d 2 2)3 + q[|(/(v2+w2) 

+ - )(r 2 l + d . 3 + d„ 3 ) 

/("VT+W 2) c 1 2 

-T^cW)]-r[j ™ 

+ d 1 3 +d 2 3 )]}, 


- (r 2 L 

/(V*+W*) 



(D.2.6) 


T = (T T T ) T 

a 3 a 3 a 3 a 3 
x y z 

where 

T = | p c , r ** [2Up - Vq - Wr], 
o a c 

x 


T a = I p C d r c 2 U[2Wx 2 + (r c 2 + 2x 2 2)q] ’ 

y 


and 

T = “ T P c , r 2 U[2Vx 0 -(r 2 + 2x 2 )r] . (D.2.7) 

a^ 4^dc Z c z 

z 

Here, X 2 is the distance from the center of the rear endcap to the 


centroid of the cylinder. 



